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Abstract 

We discuss the numerical solution of the Schrodinger equation with a time-dependent Hamilton operator using 
commutator-free time-propagators. These propagators are constructed as products of exponentials of simple weighted 
sums of the Hamilton operator Owing to their exponential form they strictly preserve the unitarity of time-propagation. 
The absence of commutators or other computationally involved operations allows for straightforward implementation 
and application also to large-scale and sparse matrix problems. We explain the derivation of commutator-free expo- 
nential time-propagators in the context of the Magnus expansion, and provide optimized propagators up to order eight. 
An extensive theoretical error analysis is presented together with practical efficiency tests for different problems. Is- 
sues of practical implementation, in particular the use of the Krylov technique for the calculation of exponentials, are 
discussed. We demonstrate for two advanced examples, the hydrogen atom in an electric field and pumped systems 
of multiple interacting two-level systems or spins that this approach enables fast and accurate computations. 
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1. Introduction 

The time-evolution of a driven quantum system is determined by the Schrodinger equation 

id,m = H{t)m (1) 

with a time-dependent Hamilton operator //(f), which one tries to solve for a given initial wave function i/'(fo) and 
times ?>?(). Prominent examples are atoms in laser fields, spins in magnetic fields or quantum dots contacted to AC 
voltage sources (see e.g. Ref. |1] for an introductory discussion). The Schrodinger equation is a special case of a 
general linear differential equation 

d,x{t) = A{t)x{t) (2) 

with time-dependent coefficients, where A(t) - -\H{t). Other examples from quantum mechanics are the Liouville- 
von-Neumann equation for the density operator p(f) or master equations for dissipative systems [2J. Analytical so- 
lutions of such equations can be found only in a very limited number of cases. In most situations one must resort to 
numerical computations. In the present paper we study an efficient numerical solution technique, which is related to 
the Magnus expansion but avoids the use of commutators. 

The propagator U(ti,t2) of Eq. ^ satisfies the initial value problem 

5„t/(fl,f2)=A(fi)f/(fi,f2), f/(fo,fo) = /, (3) 

with the identity operator, or matrix, /. The solutions x{t) of Eq. (|2]l fulfill x{ti) - U{ti,t2)x{t2)- We note the group 
property U{h,t2)U(t2,h) = U{h,h). 
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For time-independent A s A(t) the propagator is given by a (matrix) exponential 




The generalization of this expression for time-dependent A(t) is due to W. Magnus jstl- The Magnus expansion (we 
refer the reader to the recent review [Q) expresses the propagator in the form 



U{t) = U{t, 0) = exp[Q(f)] . 



(5) 



Notice that we often set the initial time to - 0. Expressions for arbitrary initial time to are obtained by the variable 
substitution 1 1-^ t + to- The operator Q(f) is given as an infinite series 



involving nested commutators of A{t) at different times. Only if [A(ti),A(t2)] - for all fi, t2, Eqs. (|5]l, (|6]l reduce to 
the simpler expression Eq. (|4). Otherwise the nested commutators provide the necessary correction terms. 

The Magnus expansion is important from a theoretical and practical point of view. In many cases the diff'erential 
equation Eq. (|2|l has an underlying Lie group structure, where the propagator U(ti,t2) is element of a Lie group 
and A(t) of the associated Lie algebra. For the Schrodinger equation, skew-hermiticity of A(t) = -iH(t) implies 
unitarity of U(ti,t2)- Violating unitarity results leads to artificial decay or growth of relevant components of the wave 
function, which spoils the stability of numerical time propagation. In particular, only unitary propagators preserve the 
normalization of the wave function. The Magnus expansion respects the Lie group structure, since the exponential 
function Q. i-> exp[Q] maps Q.{t), which as a sum of commutators of A(t) is itself a Lie algebra element, onto a Lie 
group element U(t). 

From the practical point of view, a truncation of the infinite Magnus expansion provides an approximate propagator 
U(t + 6t, f), which can be used to propagate a solution x{t) over a small time-step 6t. For an Mh-order approximation, 
where the approximation error scales as <5f'^^', all terms with or less commutators in Eq. (|6) must be kept. It 
is the virtue of the Magnus expansion that for every truncation U{t + 6t, t) is a Lie group element (whenever a Lie 
group structure is present). In this way the Magnus expansion allows for the systematic construction of geometric 
integrators ||5|-|7|, which preserve Lie group structures. 

The practical evaluation of the Magnus expansion is however rather involved. The number of terms in Q(f) is large 
already for moderate approximation order, and their calculation is complicated because of the nested commutators. 
Our starting point for better numerical algorithms are approximations of the form 



where each A, - 2„ gi^nMtn) is a (finite) linear combination of A{t) at different times f„ e [t,t + 6t] (which will later be 
chosen as Gauss-Legendre quadrature points). Such commutator-free exponential time-propagators (CFETs) preserve 
Lie group structures through the exponential form of the approximation but avoid the use of commutators. Their 
application is thus straightforward and requires only slight adjustments of existing programs for the calculation of 
matrix exponentials. No complicated scheme for the computation of nested commutators or the storage of intermediate 
results is needed. CFETs are examples for Crouch-Grossman methods ff\, and have been studied with a focus on linear 
differential equations in Refs. [8, 2]- In particular the work of Blanes and Moan [SJ, together with the review 
provided the initial motivation for the work reported here. 

In the present paper we discuss CFETs from a practitioner's point of view. Our intention is to provide a com- 
prehensive account of the theoretical background and a demonstration of the practical usefulness of this approach. A 
specific goal is the construction of optimized high-order CFETs, which can be applied to the Schrodinger equation in 
general situations where the resource consumption of naive computational approaches, e.g. a second-order approx- 
imation, would be intolerably large. To pursue these goals we first revisit the derivation of the Magnus expansion 
(Secs.|2][3]l and of the order conditions for the CFET coefficients (Sec.|4]l. A notable deviation from the Uterature is 
the replacement of a power series expansion of A(f) with an expansion in Legendre polynomials. Their orthogonality 
properties allow to simplify the presentation in two important aspects. First, the rather non-obvious fact that, effec- 
tively, only terms of order dt'^^'^ of A(f) must be taken into account for the construction of Mh-order approximations 




(6) 



Uit + St, t) = e^'e'^- 



(7) 
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is evident from the structure of the order conditions. Second, the application of Gauss-Legendre quadrature (Sec. |7) 
is straightforward, and the corresponding coefficients are obtained without additional work. We believe that our pre- 
sentation is not only simpler than others in the literature, but allows the reader to understand the derivation without 
taking unexplained aspects for granted. 

Extending previous results we construct CFETs up to order 8. Their error is analyzed theoretically in Sec. |5] 
complemented by a practical error analysis in Sec. |6] Minimization of the CFET error requires inclusion of higher 
order terms from the Magnus expansion, specifically of the N 12 + 1 -order term of A{t) for an Mh-order approximation. 
This in contrast to the error analysis for split-operator techniques found in the literature. Our improved analysis leads 
to optimized 4th- and 6th-order CFETs. Again, the use of Legendre polynomials is vital for the analysis. 

In practical applications with large Hamiltonian matrices the evaluation of the exponentials in Eq. ^ is the deter- 
mining factor for the actual efficiency. We discuss the combination of CFETs with the Krylov technique in Secs.[8j|9l 
In Sec. [To] we compare CFETs with the (f, f')-method, a Floquet-based approach. Finally, we demonstrate in Sec.fTTI 
the application of CFETs in two situations where precise results are hard to obtain otherwise, e.g. with the original 
Magnus expansion, before we conclude in Sec. [T2] The appendices give the recursion for the Magnus expansion in 
a form suitable for computer algebra computations, a short discussion of free Lie algebras and Hall bases, and the 
explicit solution of the order conditions for 6th-order CFETs. 

2. The Magnus expansion 

The Magnus expansion provides Q.{t) in Eq. (|5]) as a series 

oo 

D(f) = 2]n„(f), (8) 

71=1 

where Q„(0 is the n-fold integral of a sum of « - 1-fold nested commutators of A{f). We say that a function /(f) is 
of order t'^ if lim,^o f{t)/t^^^ - 0, i.e. its power series in t starts with . Since each integration over t increases the 
order by one, the term Q„(f) is of order t". Derivations of the Magnus expansion can be found at many places in the 
literature (cf. Ref. fT]). For our presentation, we follow Ref. (id]. The principal idea is to find an implicit equation 
relating A(t) with Q.{t), which is solved order by order for the Q„(f)- Notice that we always assume that A{t), and the 
solutions of Eq. (|2]l, are sufficiently regular to permit a local power series expansion. 

2.1. Derivation 

By definition (Eqs. (|5])), Q(f) is the solution of the implicit differential equation 

5,e"<'' = A(f)e"*'^ , 0(0) = . (9) 

To evaluate the derivative of the matrix exponential on the left hand side, consider the function f{s, t) - 5,e*"*'^. It 
fulfills the differential equation 

djis, t) = 5,5,e'"<" = 5,e'"("Q(f) = /(s, t)Q.{t) + e*"«Q(f) (10) 

with initial condition /(O, t) - 0, whose solution is given by f{s, t) = j^' e''"('> Q(f) e('-'-)"«t/r. For .? = 1, we obtain 

Wo / ,t'oWo '«! / ti^^^+^y- 

where we used the identity e^Ye^^ - 2,'^^o(l/'"')[^' ^l"' "^^^^ the iterated commutators 

[X, Y]o^Y, [X, y]i = [X, Y] , [X, y]„,+i = [X, [X, yu , (12) 

which follows, e.g., from comparison of the derivatives of 5 e'^Ye-'^ and s YZ=i)i'^" lfn\)[X, Y]„,. If the Q(f) at 
different t commute, only the first term m = in the sum contributes. Using Eq. (fTTT i in Eq. (|9]l gives 

oo 

A(f) = (5,e"») - y T-.mt), n(f)]m . (13) 
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We now insert the ansatz for the Magnus series Eq. ([8) into Eq. (fTlT l. This gives (we drop the argument f in Q(f)) 



m=0 «=1 <:=! 

(14) 



(m + 1)! , ■H' 



m-O ^ ^ ni,. ..,«,„ = ! 

To solve for Q„ we collect all terms of order f" ' . A nested commutator [Q„, , [Q,„, . . . , [Q„„,, Qj] . . . ] is of order 
rii + ■ ■ ■ + n„ + k - I in t. The only (n - l)th order term that contains 0„ is the term with m = 0. Thus, 

n-I 

Qi=A, Q„ = -y- — y [Q„,,[Q„„...,[Q„,„,Q,]...] . (15) 

(m + 1)! , 

m=l n],...,n,if,k>l 
niH h«,„+/:=n 

Notice that all sums are finite (the last term in the sum over m, for m = n - 1, is (l/n!)[Qi, Qi]„_i). A final integration 
gives the explicit expressions 

Ait')dt' , a,(o = -27 — ?TT Z Uf'["«xo,[t^«.(f'),.-.,[a,„(f'),fi^(0].-.] , (16) 

Q m=I «i,...,n,„,/t>l 

n 1 H — hn,„ +/:— » 

which allow for the recursive calculation of the Q.„{t). As stated before, every term in Q„(f) involves an n-fold integral 
of an n - 1-fold nested commutator of A(t). We obtain explicitly, up to order f^, 

Q(f) = (f) + ^2(0 + Q3(f) + Oit"^) , with 



Qi(f)= r dhAiti), Q2(0 = ^ r dt2lAiti),Ait2)], 
Jo ^ Jo Jo 

Q3(0 = ;^ r fl'fi f'fl'fi r't!'f3[A(fi),[A(f2),A(f3)]] + [[A(f,),A(f2)],A(f3)] . 
o Jo Jo Jo 



(17) 



It is convenient to write the Q„(0 as time-ordered integrals [10]. This requires additional manipulation of the 
integration domains of the terms found by straightforward integration in Eq. ( fTSI l. It is possible to derive a systematic 
recursion (see Appendix A , which is very useful for symbolic calculations on a computer 

An alternative route to solve Eq. (fT3T l is to note that the commutator expression on the right hand side involves the 
Taylor expansion of the function (e' - l)/x = Tjm=o x'"/iin +1)1. Solving for Q(f) is thus possible using the inverse 
function, where the Bernoulli numbers B„ appear as the Taylor coefficients in x/(e* - 1) = ET^o ^"fj- After a few 
additional manipulations one obtains again a recursive definition of the Q.„(t) (see e.g. Ref j^). Our experience is 
that the present approach is better suited for an algorithmic implementation. Interestingly enough, it avoids the use of 
Bernoulli numbers. 



2.2. The Baker-Campbell-Hausdorff formula 

A special case of the Magnus expansion is the Baker-Campbell-Hausdorff (BCH) formula 

e^e^ = exp[x + y + ^[X, Y] + ^[X, [X, Y]] - ^[Y, [X, Y]] - ^[F, [X, [X, Y]]] + ...]. (18) 

We note that the left hand side of this equation is the exact propagator for a stepwise constant A(f), with A(f) - Y for 
< f < 1, A(f) - X for I < t < 2. Inserting this A(f) into the recursion Eq. (fTST i for D.(t) provides the exponential on 
the right hand side. In a similar spirit, we can obtain the BCH formula for several exponentials 

" 1 

-exp[2x, + - ^ + . . .] . (19) 

/=1 l<i<j<s 
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3. Approximate Magnus propagators 

By construction, the Magnus expansion is an expansion in orders of t. It thus provides a systematic way to obtain 
Mh-order approximations U'-^HO - exp[2^=i ^ri(t)]^ which coincide with the exact propagator U(t) for all terms of 
order or less, from direct truncation of the infinite series Eq. (|8). Notice that we call a function /(f) an Mh-order 
approximation of another function g(t) if the difference f(t) - g(t) is of order t^^^ . 

The expression for Q„(f), given through Eq. ( fTSI l. involves «-fold integrals. These can be simplified since each 
integral needs to be evaluated only up to order f'^^' for Mh-order approximations. Starting from an expansion of A(t) 
in orders of f, all multi-dimensional integrals in Eq. (fTTT i can be replaced by one-dimensional integrals. In the literature, 
it is common to expand A(f) in powers of f (or centered powers (f - 6t/2)" for a given time-step 6t). Contrary to these 
treatments, we use an expansion in Legendre polynomials. Although both expansion are principally equivalent, the 
choice of Legendre polynomials proves itself useful because of their orthogonaUty properties. 

3.1. Legendre expansion ofA(t) 

The (shifted) Legendre polynomials P„{x) are defined for « = 0, 1, 2, . . . through the recurrence 

Po(x) = 1 ,Pi(x) = 2x - 1 ,P„+i(x) = ^^(2x - l)P„(x) - -^P„_,{x) . (20) 

n+\ n+\ 

By definition, P„{x) is a polynomial of degree n. Explicitly, 

Piix) ^6x^ -6x+l , P2(x) = 20x^ - 30x^ -H 12jc - 1 , P4(x) = 70/ - 140x^ + 90x^ - 20x + 1 . (21) 



The polynomials P„{x) are symmetric with respect to x = 1 /2, i.e. 

P„(l-x)^(-irP„(x) . (22) 

Furthermore, they form a complete set of orthogonal functions on the interval [0, 1], with scalar product 

1 

P„rix)P„(x) dx = -6,n„ . (23) 

J 2« H- 1 



In particular, p(x)Pn(x) dx -0 for every polynomial p(x) of degree less than n. 

We now fix a time-step 6t, for which an approximate Mh-order propagator U^'^HSt) = U^'^HSt, 0) should be 
constructed. The function A( f) is expanded on the interval [0, 6t\ in a series of Legendre polynomials 

1 ^ 

A{t) ^ -Yj^nPn-i{jy 0(,61^^') {Q<t<5t). (24) 

;7=1 

Notice the index shift of P„_i versus A„. The (matrix-valued) coefficients are obtained as 

St 1 

A„ = (2n-1) ^ A{f)P„^x[j^dt^{2n-\)5t J A(x6t)P„^i(x)dx . (25) 


To see that A„ is a term of order 6f, compare this expansion with an expansion A(f) = I],„>i fl/ji?'""' in pow- 
ers of t. Since Pnix) is orthogonal to all x™ with m < n, we see from Eq. (l25T l that A„ starts with the term 
6t a„(xSt)"^^ Pn^i(x) dx of order St". In particular, it is a linear combination only of with m >n. 
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3.2. Legendre expansion of Q.(t) 

If we insert the expansion Eq. ( l24l i of A(f) into the recursion Eq. (fTSI l we obtain Q.(6t) as a sum of nested commuta- 
tors of the expansion coefficients A„. A nested commutator [A„,, ... ,A„_J is oforderni H i-n„, in 6t. Theprefactor 

of this term is obtained as the n-fold integral ^(mi, . . . , «„,) = Jq dx\ . . . j^^'"'' dx„P„^-l(xl) ■ ■ ■ Pn,„-\{xm), which is a 
rational number independent of A(f). For example, 

Q2(ft) = -^ r dh r'flff2[y A„,p„,_i(^),y A„,p„,_i(|)] 

" ^ y! I r "^^1 r '^-*^2^'n,-i(-*:i)^'n2-ite)| [A„,,A„J (26) 



,"2 2 

-1 



a!;c2^o(^i)Afe)-^i(^i)i'ofe)| [Ai.Aj] = -^[Ai,A2] 



Notice that the only non-zero contributions in the second line come from n\ - \,n2 - and ni - 0,n2 - 1, since the 
integral of P„, (xi)P„, (jC2) vanishes in all other cases. This hints at a general pattern to be discussed below. 
Collecting all terms up to order 6t'^, we find 

Q(50=Ai-i[Ai,A2] 
o 

+ -^[Ai, [Ai, A3]] - ^[A2, [Ai,A2]] + -^[Ai, [Ai, [Ai, A2]]] - 4t[^2, A3] 

- ^[Ai, [Ai, [Ai,A4]]] - ^[[Ai,A2], [Ai,A3]] + ^[A2, [A,, [AuA^]]] - ^[Aj, [A2, [Ai,A2]]] 
^ -[A3, [Ai, [Ai, A2]]] - :r^[Ai, [Ai, [Ai, [Ai,A3]]]] - J-[[Ai, A2], [Ai, [Ai, A2]]] 



2520 2520 7560 
1 1 



[A2, [Ai, [Ai, [Ai,A2]]]] - TT^[Ai, [Ai, [Ai, [Ai, [Ai,A2]]]]] + 0(5?'°) . 



2520 15120 

The first line contains the 4th-order terms, the second line the 6th-order terms, and the remaining lines the 8th-order 
terms. This expression for Q((5f) avoids multi-dimensional integrals. 

5.5. Properties of the expansion 

As seen above for D.2{5t), only few out of the many possible commutators contribute to Q.{5t). In Eq. ( l27l l several 
nested commutators of order 5t'^ or less are missing, e.g. the terms A2, . . . , Ag or [Ai, A3], . . . , [Ai, Ag]. This is a 
consequence of two general properties of the expansion that result in a zero prefactor ^(ni, . . . , n,,,) of [A„,, . . . , A„,J. 

(PI) Time-reversal symmetry U{6t, 0) - U(0, 50 ' of the propagator implies that Q(f) changes sign if A(f) is replaced 
with -A{6t - t). According to the parity Eq. (l22t of the Legendre polynomials it follows that even order terms in the 
expansion, i.e. terms [A„j, . . . , A„,J with even n\ + ■ ■ ■ + n„„ vanish. This follows also from the calculation of the 
prefactor ^(ni, . . . , «,„) as an n-fold integral: Each of the inner integrations over X2,...,Xn changes the parity of the 
integrand. The parity also changes by multiplication with a polynomial P„ for odd n. Hence, the integrand in the final 
integration over x\ has odd parity for odd {ni - Y) + ■ ■ ■ + - 1) H- (-1)™"', i.e. if «i + ■ ■ ■ + n^ is even. Then, the 
integration gives zero and the respective term vanishes in Eq. ( IZTl i. 

ff 2j As a consequence of the orthogonality of the Legendre polynomials a term [A„j , . . . , A„ J vanishes if some index 
nii exceeds the sum of the others by two, i.e. n^ > 1 + f^i" a = 1, . . .,m. To see this change the integration 

order in the integral for the prefactor ^(ni , . . . , «„,) such that the outermost integration is over x^. This final integration 
is of the form dx^P „^-\{xk)p{xk), where the polynomial p{xk) results from the previous m - 1 integrations of the 
other polynomials P„__i(x,). The degree of p{xk) is at most 2,^^ n,. If, by assumption, this sum is smaller than nk - 1 
the final integral is zero since is orthogonal to polynomials with smaller degree. Now suppose [A„,, . . . , A„„J is 
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Ai, A3, [Ai,A2], [Ai,A4], [A2,A3], [A3,A4] 
[Ai,[Ai,A3]], [A2,[Ai,A2]], [A2,[Ai,A4]], [A2, [A2, A3]] , [A3, [A,, A3]] , [A4,[Ai,A2]] 
[Ai,[Ai,[Ai,A2]]], [Ai,[Ai,[Ai,A4]]], [A2, [Aj, [Aj, A3]]] 
[A2,[A2,[Ai,A2]]], [A3,[Ai,[Ai,A2]]], [[Aj, A2], [Aj, A3]] 
[Ai, [Ai, [Ai, [Ai, A3]]]] , [A2, [Ai, [Ai, [Ai, A2]]]] , [[Aj, A2], [Aj, [Ai, A2]]] 
[Ai,[Ai,[Ai,[Ai,[Ai,A2]]]]] 

Table 1: The 22 odd elements up to order 8 of the Hall basis with generators Ai, . . . ,A4. 

a term of order A^, and one index rii^ > N/2. Then, > 2,- - rik + Yjo^k > ■^/2 + Tjo^k "i' or Yjo^k < N/2. The 
above condition appHes, and it follows that this term gives no contribution. 

Both properties considerably simplify the derivation of approximate propagators since they reduce the number of terms 
in the expansion Eq. dZTl l. According to (PI), approximately only half of the commutators contribute. In particular, 
an expansion including all terms up to some odd order 6t^ is automatically correct up to order df'^^'. According to 
(P2) there are no contributions from terms such as A,„ for m > 2 or [A,„, A„] for \m - «| it 1. It also explains why in the 
Mh-order expansion of 0(f) only commutators of terms Ai, . . . ,Aff/2 occur (i.e. Ai, . . . ,A4 in Eq. (|27]|). This has the 
remarkable consequence that for the construction of Mh-order propagators terms A„ for n > N/2 + 1 can be neglected 
even \fn<N. Notice that property (PI) is shared by an expansion in centered powers (f - 6t/2)", while (P2) requires 
orthogonality of the Legendre polynomials. This fact motivated our use of a Legendre expansion of A(f) instead of 
the apparently simpler Taylor expansion. 

3.4. Uniqueness of the expansion: Hall basis 

The expression for Q.{6f) in Eq. (|27] | is not unique. Non-trivial identities between nested commutators, e.g. the 
Jacobi identity [A, [B,C]]h-[B, [C,A]]h-[C, [A,B]] = 0, allow to replace one commutator by others. To compare nested 
commutator expressions by equating the coefficients we must therefore first eliminate the ensuing linear dependencies. 
Technically, this amounts to calculations using a vector space basis of the free Lie algebra generated by the A„. Since 
every nested commutator is a unique linear combination of the basis elements, uniqueness of the entire Magnus 
expansion is achieved. 

A systematic construction of free Lie algebra bases is provided by a Hall basis |11]. Algorithms exist for the 
rewriting of nested commutators in terms of the Hall basis elements, and for their enumeration. The number of Hall 
basis elements grows rapidly with the maximal order considered. As listed in the following table, 

order N 2 4 6 8 10 

full set of elements 2 7 22 70 225 

relevant according to (PI), (P2) 1 2 7 22 73 

there are 70 elements up to order 8 in the Hall basis. As a consequence of the two properties (PI), (P2) from Sec. 13.31 
only the 22 elements in Table[T]are relevant for our purposes. Notice that in Eq. dZTl l the elements A3 and [A 1 , A4] from 
the Hall basis are missing according to (P2), but yield order conditions for the CFETs as discussed in Sec. 14. II We do 
all calculations using the 22 Hall basis elements, rewriting commutators as necessary, e.g. [Ai, [A2, {A\, [Ai,A2]]]] — 
[A2, [Ai, [Ai, [Ai, A2]]]] + [[Ai, A2], [Ai, [Ai, A2]]]. 



4. Commutator-free exponential time-propagators 

The expansion Eq. (l27T i of Q.(6t) still contains nested commutators. An Mh-order commutator-free exponential 
time-propagator (CFET) is based on the ansatz 

f7[,^'(5f) = ^^"'e"---^^"' , (28) 
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where each of the s exponentials Q, is a linear combination 

N 

Q, = 2/;-,„A„ (29) 

n=l 

of the Ai, . . . ,An from the Legendre expansion Eq. (l24l l of A(t). The CFET is completely determined through the 
coefficients /j „, which are fixed once and independently of the concrete A(t) used in a calculation. The practical 
evaluation of Eq. (l28T l. avoiding commutators and multi-dimensional integrals, is considerably simpler than for the 
original Magnus expansion. It will be discussed in more detail in Sec. [8] 

Effectively, Eq. (l28T l is the exact propagator for an auxiliary problem with a fictitious, stepwise constant A(f). The 
CFET coefficients / „ must be determined in such a way that the replacement of the complicated time-dependent 
problem by the simpler auxiliary problem introduces only an error oc 5?^^+', independently of A(f)- Now consider an 
A(f) - x(t)X + yit)Y which is the sum of two contributions X, Y. This situation arises, e.g., for a particle moving in 
a time-dependent field. By construction, each Q, - XiX + yiY itself is a sum of X, Y, with constant Xi,yi replacing 
x(t),y(t). Therefore, the CFET describes again a particle moving in a field, and thus preserves the principal physical 
situation. Notice, however, that fictitious negative time-steps can occur The analogous statement does not hold for 
the original Magnus expansion involving commutators of X, Y. 

The simplest example of a CFET is the 2nd-order midpoint rule 

^CF2:i(^') = exp[Ai] = exp [ J dtAit)] ^ exp[(5fA(5f/2)] , (30) 

corresponding to i = 1 and fu - 1. The second exponential is identical to the first according to the definition 
Eq. ( |25] ) of Ai. The last exponential is obtained by approximation of the integral through Gauss-Legendre quadrature 
(addressed later in Sec. |7]i, which here reduces to evaluation of A(t) at the midpoint f = 6t/2. 



4.1. Derivation of order conditions 

The construction of higher-order CFETs is substantially more difficult, and a systematic procedure is missing. We 
adopt the following strategy: Starting from the CFET ansatz Eq. ( l28T l. the BCH formula ( fT9l l allows us to combine the 
s exponentials until we obtain U^^^{5t) - with 

/\ \^ n A V fi.\fi,2~ fi.yfi,2 , . . , 

" = Zj ~2 [M,A2\ + ... (31) 

1=1 \<i<j<s 

The Cl has to be compared with Q.{St) from the Magnus expansion Eq. (l27l l. demanding equality of terms of order ^f^^ 
or less. Working in a Hall basis, this implies equality of their prefactors which results in equations for the coefficients 
/i,„, the so-called order conditions. Specifically, we find 

s 

2]/;> = 5«,i (32) 

1=1 

arising from the terms A„, and from [Ai, A2] 

X klh2-fiAk2--\- (33) 

l<!<;<i 

For higher-order commutators, the derivations become increasingly cumbersome, and calculations are best delegated 
to a computer Since standard computer algebra systems are less useful for calculations in non-commutative algebras 



we used self- written programs that perform the Lie algebra manipulations, based on algorithms from Ref. H12 I. 

Counting all Mh-order elements in the Hall basis (Sec. 13.4b . we see that the number of order conditions is 22 (70) 
for order 6 (order 8), and thus appears to be too large for a practical solution of the multivariate polynomial equations 



that arise. As we found in Sec. l3.2l several commutators do not appear in Q.{6t) as a consequence of the two properties 
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4th-order 


2 exponentials 


3 exponentials 


CF4:2 


CF4:3 


/i,i = l/2 /i,2 = l/3 /2,i=0 


/u = 11/40 /i,2= 20/87 /2,i=9/20 



Table 2: Coefficients for 4ffi-order CFETs wiffi 2 and 3 exponentials. Notice that CF4:3 is not recommended for use (cf. Sec. 16.31 . 

(PI), (P2). The key observation is that the corresponding order conditions can be satisfied by a suitably restricted 
choice of the according to the following two rules. 

(Rl) Since passing from U^^p{5t,U) to t/*^'(0, 5f) ' changes the sign of A„ by (-1)", a CFET complies with time- 
reversal symmetry if the coefficients obey 

f.-M,n = (-ir^'/i> . (34) 

For a time-symmetric CFET it thus suflices to specify the // „ for / < i/2, i.e. for the first half of the exponentials 
e"', and choose the remaining coefficients according to Eq. ( l34l i. For odd s, the coefficients /(s+i)/2,« of the central 
exponential / = {s + l)/2 must be specified for odd n only, while they are zero for even n. With this constraint the 
order conditions for even order terms, which do not contribute to Q(5f) according to (PI), are automatically satisfied. 
(R2) Property (P2) states that up to order only terms A„ with n < Nj2 contribute to f2(5f)- The order conditions 
involving higher-order A„ can be satisfied simply by setting - for n > N/2: Since all coefficients are zero the 
corresponding commutators drop out entirely. The remarkable implication is that an Mh-order CFET can be built 
already from the terms Ai, . . . ,An/2- We note that this property is intrinsically connected with Gaussian quadrature 
using orthogonal polynomials (cf. Sec.|7]i. It becomes obvious working with Legendre polynomials, while it requires 
sophisticated additional arguments in general fl. 

By rule (Rl) the number of relevant coefficients and order conditions is reduced approximately by one half. For this 
reason we consider only time-symmetric CFETs. Notice that a symmetric Mh-order CFET is automatically of order 
-I- 1, if is odd. Rule (R2) implies that the summation index n in Eq. ( |29] l only has to run from 1 to N/2. We will 
later relax this rule to allow for minimization of the error, which requires inclusion of the term Aai/2+i . 

With both rules, the number of order conditions is significantly reduced, to 2, 7, 22 for = 4, 6, 8 CFETs (cf. 
the Table in Sec. 13.41 ). On the other hand, a symmetric A^th-order CFET with s exponentials has [sN/A] coefficients 
(rounding down to an integer). The counting shows that 5 exponentials (11 exponentials) are needed for a 6th-order 
(8th-order) CFET. Only in exceptional cases solutions with less exponentials exist, e.g. CF6:4 in Table[3] 

4.2. Fourth-order CFETs 

We consider 4th-order propagators (A^ = 4) with three exponentials (s - 3), of the form 

U^^^(6t) = exp[/i,iAi + /i,2A2] exp[/2,iAi] exp[/i,iAi - /i,2A2] . (35) 

As explained before (cf. Eqs. (l32l l. (l33Tl). we get the two order conditions 

1 =2/i,i +/2,i , 

1 (36) 

~7 - ~(/l,l + /2,l)/l,2 • 
O 

The first arises from the term Ai, and the second from the term [Ai, A2]. In accordance with the above counting of 
terms, we have 3 coefficients and 2 order conditions. Using /2j as the free parameter, we find 

1-/21 1 

Corresponding coefficients are listed in Table |2] The parameter /2 J will later allow for optimization of the propagator 
(see Sec. 15.2b . Setting /2,i = 0, we obtain the unique 4th-order CFET with s - 2 exponentials 

U§,.^(6t) = exp [^Ai + JA2] exp [^Ai - JA2] . (38) 

The notation used here and in the following is CFN:s for an A^th-order CFET with s exponentials. 
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6th-order 



4 exponentials 
CF6:4 

_ 1 (5400- 600 V6)'^^ (^(9 + V6))'^' 
/I.I - 2 60 ^2- 32/3 ^ 

/i,2 = /i,i - /i,3 = Yo^Yo/l 

/2,l = i-/l,l /2.2 = i(l- 4/1,1+ 2/2j) /2,3 = -/l,3 



/3,I = 


/3,2 - 





/3,3 - 





5 exponentials 


CF6:5 


/i.i = 0.16 


/l.2 = 


0.14587456942714338561 


/l,3 = 


0.11762370828143015682 


/2 1 = 0.38752405202531 186588 


/2,2 = 


0.15089113704380764664 


/2,3 = 


-0.12805075909013044594 


/3,i= 1-2/2,1-2/1,1 


/3,2 = 





/3,3 = 


-2/2,3-2/1,3 


CF6:5b (cf. Ref. [8]) 


/i,i= 0.2 


/l,2 = 


0.1746879190177786220 


/l,3 = 


0.12406375705333586606 


/2i= 0.34815492558797391479 


/2,2 = 


0.1068765450953683 


/2,3 = 


-0. 1 3902 1 3 1 3323765096675 


/3,1= 1-2/2,1-2/1,1 


/32 = 





/33 - 


-2/2,3 - 2/1,3 


6 exponentials 


CF6:6 


/u = 0.16 


/l.2 = 


0.15101538937746543493 


/l,3 = 


0.13304616813239630479 


/2,i = -0.22738164742696330169 


/2,2 = 


-0.087654259755 1 1543 1662 


/2,3 = 


0.069919836812656575583 


/3,i= l/2-./i,i-/2,i 


/3,2 = 


0.21035154512209824847 


/3,3 = 


-/l,3 -/2,3 



Table 3: Coefficients for unoptimized 6th-order CFETs wiffi s = 4,5,6 exponentials (coefficients for optimized 6th-order CFETs are given in 
Table|6j. The CFET CF6:5b corresponds to the coefficients of the propagator i/ij^' from Ref. Q]. 



4.3. Sixth-order CFETs 

For 6th-order (A^ - 6), we consider propagators with s - 6 exponentials. The 9 coefficients /■,„, for 1 < /, n < 3, 
must satisfy 7 order conditions corresponding to the 7 Hall basis elements 

Ai , A3 , [Ai, A2] , [A2, A3] , [Ai, [Ai, A3]] , [A2, [Ai, A2]] , [Ai, [Ai, [Ai.Aj]]] (39) 

from Table[T] We note that two coefficients can be chosen as a free parameter 

An explicit solution of the order conditions is possible to a large degree, and simple explicit expressions for the 



coefficients can be obtained in some cases (cf. Appendix C 1. Setting /3,2 = 0, the two central exponentials can be 



combined, resulting in propagators with 5 = 5 exponentials and a single free parameter Surprisingly, there is also a 
solution with /3,i - f2^2 - /3,3 = 0, giving a 6th-order CFET with only 4 exponentials (see CF6:4 in Table|3ll, although 
there are less coefficients than order conditions. We do not know whether the existence of this solution is accidental, 
or hints at a general redundancy pattern of the equations. For practical purposes, the CFET CF6:5 from Table |3]is 
most relevant, since it has small approximation error Further optimized 6th-order CFETs will be obtained in Sec. 15.31 

4.4. Eighth-order CFETs 

The 22 order conditions of 8th-order CFETs correspond to the entire set of commutators from Table [1] Exactly 22 
coefficients exist for i = 1 1 exponentials. Due to their complexity, the order conditions can only be solved numerically. 
Several solutions were computed using a root finder based on the Newton iteration (v^. Severe ill-conditioning of the 
equations required the use of high-precision arithmetics, based on the MPFUN package 1IT4I1 . and repeated restarting 
of the Newton iteration. The coefficients of an 8th-order CFET with small approximation error, selected from about 
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8th-order: 1 1 exponentials 



CF8:11 



f, , 






/1, 2 - 




/l,3 




0.119167378745981369601216 


A4- 


0.0686 19226448029559 107538 


/2.I 




0.3794208075 1600543 1504230 


A2 - 


0. 148839980923 180990943008 


/2,3 




-0. 1 15880829 1 8662807502 1088 


flA - 


-0.188555246668412628269760 


/3.I 




0.4694593066440505730 1 7994 


hi- 


-0.379844237839363505173921 


/3,3 




0.0228988 14729462898505 141 


A4 = 


0.571855043580130805495594 


Ai 




-0.448225927391070886302766 


A2 - 


0.362889857410989942809900 


A3 




-0.022565582830528472333301 


f4A - 


-0.544507517141613383517695 


Ai 




-0.293924473106317605373923 


/5.2 = 


-0.026255628265819381983204 


/5.3 




0.09676 1 509 1 3 1620390 1 00068 


/5.4 = 


0.000018330145571671744069 


Ai 




0.447 1 095 1 05867986 1 4 1 20629 


A3 = 


-0.200762581179816221704073 



Table 4: Coefficients for an 8th-order CFET wiffi 1 1 exponentials. 



50 computed solutions of the order conditions, are given in Table |4] A systematic search of the coefficient space was 
not possible. 

5. Theoretical error analysis 

The CFET error is determined by the difference = Q - Q between the exact Q(5f) from the Magnus expansion 
Eq. (l27l l and the approximate O from Eq. (ISTT i. The theoretical error analysis aims at minimization of the error term 
in the general situation, where no specific information about A(f) is available. 

5.7. General considerations 

By construction the error term is of the form x = HkiPk - Ck)Ck, where the Ca are the + 1 -order commutators 
from the Hall basis, the are polynomials in the coefficients /j „ such as in Eq. (ISTT l. and the q ffie constant prefactors 
from Eq. dZTl l. The size ofx can be measured with a matrix norm || ■ ||. It is 

IWI = II J]^Pk - Q)Q|| <Yj\Pk- ■ WCkW ■ (40) 

k k 

In concrete situations, \\x\\ depends not only on the size ||Q|| but also on the amount of dependency between 
diff'erent Ck, which is responsible for the diff'erence between left hand and right hand side of the above inequality. In the 
general case we may not assume ffiat the difference is small. Accidental cancellations, i.e. \\{pk-Ck)Ck+ipi-ci)Ci\\ ~ 
for ?ik + I, are typical. Optimization of a CFET, that is minimization of \\x\\ through variation of the coefficients /j,,, 
thus requires that all \pk - Ck\ become simultaneously small. Only then, the error can be expected to be small in 
the general case. Such optimized CFETs are universally applicable and perform equally well in different situations. 
Optimization will be achieved for 4th- and 6th-order CFETs, listed below in Tables |5] |6] For 8th-order CFETs, 
optimization is not practicable due to the complexity of the order conditions. 

An important point, which seems to have been missed in the literature, complicates the error analysis in compar- 
ison to split-operator techniques. While rule (R2) in Sec. 14.1 [ states that the terms A„ for n > N/2 can be disregarded 
in the construction of an Mh-order CFET, the error term x contains a contribution from Aiv/2+1 since the prefactors 
Ck of the corresponding commutators are non-zero in Eq. (l27i . For 4th-order, this applies to the terms [A], [Ai, A3]] 
and [A2, A3] involving A3. In contrast to the basic construction of higher-order CFETs with (R2), CFET optimization 
requires explicit inclusion of An/2+i- Therefore, the optimized 4th- and 6th-order CFETs include non-zero coefficients 
for the A3 or A4 term, respectively. Additional order conditions, e.g. Yj']=i .fi,N/2+\ = arising from the An/2+\ term 
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4th-order (optimized) 



3 exponentials 

CF4:30pt ~ 
^ 11/40 /i,2 ^ 20/87 /i.3 ^ 7/50 Ai ^ 9/20 /2,3 ^ -7/25 



Table 5: Coefficients for the optimized 4th-order CFET CF4:30pt with 3 exponentials (Eq. )43t ). The unoptimized CFET CF4:3 from Table|2]is 
obtained by dropping the A3 terms from CF4:30pt. 



itself, must be accounted for. We note that for split-operator techniques fT?], where essentially the full Magnus prop- 
agator e"^'^ is replaced by the term e(^+^)' for a time-independent A{t) = X + Y, the equivalent coefficients = 0, and 
no additional provisions are necessary. 

Since the error term x is of order iSf'^^' we must ask whether also terms A„ for N/2 + I < n < N + I need be 
considered. Property (PI) states that nested commutators involving these terms do not occur in Eq. (l27t up to order 
-I- 1, i.e. the corresponding prefactor c^^ = in Eq. ( l40t . Since we have set the coefficients of these A„ to zero by 
rule (R2), they do not contribute to x and need not be considered. Notice again that the use of Legendre polynomials 
in (I24I 1 simplifies the derivation: With a power series expansion all terms up to order 6t^'^^ would explicitly contribute 
tox, and minimization of \x\ would result in a number of additional though redundant equations. 

5.2. Optimized fourth-order CFETs 

For a 4th-order CFET including the A3 term, we make the ansatz 

U^ck'^dt) = exp[/i,iAi + fi^zAi + /i,3^3] exp[/2,iAi + fi^^A^] exp[/i,iAi - /i,2A2 + /i,3A3] , (41) 

in extension of Eq. ( l35T l. The previous order conditions still apply, and /u , /i_2 are given by Eq. ( l37T i. The new order 
condition arising from the A3 term is 2/13 -1-/2.3 = 0, which gives one additional free parameter /2,3 with /13 = -/2,3 /2. 
With these choices, we obtain for the error term 



/I 1 + 2 /2 I \ /I /2 1 \ 

y = - O =( )[A2, [Ai, A2]] + { =^ )[Ai, [Ai, [Ai, A2]]] 

^ ^60 54(1 + foi)^' M440 288'' 



54(l+/2,i)- 



(^(1 - . ^)[^i>[^i'^3]] .(g^ . ^)[A2,A3] .0(..^) 



(42) 



It has four contributions corresponding to the 5th-order exponentials in the second line of Eq. ( l27l i. 

In Fig.[T](left panel) we show exemplarily \pk - q| for k = [A2, [Ai, A2]] and k - [Ai, [A], [A], A2]]] as a function 
of /2_i. The optimal choice is close to /2,i ~ 0.45 = 9/20. If we only try to minimize the contribution from these two 
commutators, neglecting the A3 terms, we thus obtain the CFET CF4:3 from Table |2] We will see below in Sec. 16.31 
that this CFET is far from being optimal. Full minimization of x through variation of /2,i and /23, including the A3 
terms, results in /2_i * 0.45, /2,3 ~ -0.28. For an optimized 4th-order CFET we thus propose the choice /2_i = 9/20, 
/23 = -7/25, which results in (cf. Table|5]l 

1 1 20 7 9 7 11 20 7 

t/cF4:3Opt(50 = exp [-A, + -A2 + -A3] exp [-A, - -A3] exp [-A, - -A2 + -A3] . (43) 



5.3. Optimized sixth-order CFETs 

Extending 6th-order CFETs with 6 exponentials (Sec. 14.3) by inclusion of the term A4 provides us with three 
additional coefficients /i,4, /2,4, 73,4. The new order condition 

= /i,i/i,4 + 2/2,1/1,4 + ./2,I./2,4 + 2/3,1/1,4 + 2/3,1/2,4 + ./3,l/3,4 (44) 

arising from the commutator [Ai, A4] fixes the value of /3.4. Notice that the term A4 itself does not lead to a new 
order condition, since it is even and the associated Cjt = = by rule (Rl). Using the explicit solution of the order 
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Figure 1: Left panel: Contributions \pi, - q| to the 5tli-order error term x of 4tli-order CFETs (Eq. J42t ) as a function of tlie free parameter /^ i, 
for commutators [A2, [Ai, A2]] (black) and [Ai, [Ai, [Ai, A2]]] (red). Both contributions become small in the vicinity of /2.1 = 0.45. Right panel: 
Maximum ;^'123 of the contributions - q| to the 7th-order error term x for 6th-order CFETs, excluding contributions from the A4 term. As 
indicated in the figure, the maximal contribution comes from either the [[Ai, A2], [Ai, A3]] (for /i 1 < 0.16) or the [A3, [Ai, A3]] term (fi i > 0.16) 
The upper panel indicates the number of solutions of the s = 5, N = 6 order conditions. 



conditions (cf. Appendix C 1, we can minimize the error term x through variation of the four free parameters fi i, 
./3,2, ./i,4, fiA- For i = 5 exponentials, we set /3_2 = /3,4 - 0. To illustrate the typical behaviour, we show in Fig. [1] 
(right panel) the partial error ;i'i23 - max Ipi^ - q| including only the contributions from commutators Ck without the 
A4 term. It depends on the single parameter f[ [. Optimal choices occur around /i 1 x 0.16, corresponding to the 
CFET CF6:5 from Table|3] This also provides partial justification for the CFET CF6:5b with /u = 0.2 from Ref. 
Inclusion of the A4 term and subsequent minimization of the associated error contribution, keeping /u fixed, results 
in the improved CFET CF6:5Imp. The full minimization of \x\ with free variation of all parameters results in the 
optimized 6th-order CFETs CF6:50pt and CF6:60pt listed in Tabled 



6. Practical error analysis 

The theoretical error analysis results in optimized CFETs, whose error is expected to be small in the general case. 
In a concrete situation dependencies between the nested commutators in the error term may lead to different results. 
To confirm the validity of the theoretical error analysis we study the CFET error for a driven two-level system. Further 
issues of practical relevance concern the choice between CFETs of different order, and the time-step selection. 

6.1. Time-stepping and effective error 

In the standard time-stepping approach, the approximate propagator U{t) over longer propagation times is con- 
structed as a product of short-time CFETs u''^p(t+6t, t). Equivalently, a concrete solution x{t) is repeatedly propagated 
over a small time step 6t. The propagator f7(r) for the maximal propagation time T is a product of A^, - T/St CFETs. 
Intermediate results are obtained at multiples of St. 

The accuracy of time-stepping is controlled through the size of 6t. For Mh-order CFETs f7<,p'(f + 6t,t), the 
error contributed by each one scales as 5t^*^ . Due to accumulation of errors, the propagation error after A^, steps is 
e = cNsSt^^^ - cTSt^ with an error constant c which depends on the concrete situation. To achieve a given accuracy 
requires a time step St < {e/cTY^^ for a maximal acceptable error e. Usually, St « T. The computational effort 
is dominated by the A^j-fold evaluation of the s exponentials in Eq. (|28] |. It is thus proportional to sNs = sT/St > 
^j\+\IN ^-\IN j.jjg effective error constant 

c = sc'l"" . (45) 

This quantity determines the efficiency of time-propagation with an Mh-order CFET with s stages. As a rule of thumb 
we note the relation 

effort oc error'^'^ x time . (46) 
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6th-order (optimized) 



5 exponentials 



CF6:5Imp 



/i,i = 0.16 


hi 




0.14587456942714338561 


/l,3 




0.11762370828143015682 


/2i= 0.38752405202531186588 


fi,i 




0.15089113704380764664 


/2,3 




-0.12805075909013044594 


/3,1= 1 - 2/2,1 - 2/u 


h,2 







/3,3 




-2/2,3-2/1,3 


/i,4= 0.074 


flA 




-0.212530296697694739551 


/3,4 







CF6:50pt 


/i,i = 0.1714 


fl.2 




0.15409059414309687213 


/l,3 




0.11947178242929061641 


/2,i = 0.37496374319946236513 


fl.l 




0.13813675394387646682 


/2,3 




-0.1 3090674649282935743 


/3,1= 1-2/2,1-2/1,1 


'fx2 







/3,3 




-2/2,3 - 2/1,3 


/i,4 = 0.07195 


/2,4 




-0.21123356253315514306 


/3,4 







6 exponentials 


CF6:60pt 


/i,i = 0.3952 


/1.2 




0.35629343479227292880 


/l,3 




0.27848030437681878641 


/2,i = -0.22432144875476807927 


/2,2 




-0.19935407393749030416 


/2,3 




-0.15625650102884866893 


/3,i= 1/2-/1,1-/2,1 


/3,2 




0.1145 


/3,3 




-/l,3 -/2,3 


/i,4= 0.1579 


/2,4 




-0.09512 


/3,4 




-0.16475168057141371958 



Table 6: Coefficients for optimized 6th-order CFETs wiffi .s = 5, 6 exponentials. In each case, the last row gives ffie coefficients for the A4 term. 
The CFET CF6:5Imp is obtained from CF6:5 (Table|3) thi'ough separate minimization of the A4 error contiibutions. 



6.2. Driven two-level system 

Our test problem is a driven two-level system, realized, e.g., by a spin 1/2 in a magnetic field B{t) - (Bjc{t), Byit), B,{t)). 
In the eigenbasis of the z-component of angular momentum, the Hamilton operator is given by the matrix 



mt) = 2 Z ^^^^^'^^ = 2 



k=x,y,z 



B,(t) - iBy(t) 



2 \BAt) + iBv(f) -B,it) 



with the standard Pauli matrices III6II 



cr , = 



1 

1 



-i 

1 



1 
-1 



(47) 



(48) 



For particular choices of B(t) the propagator can be expressed in simple, closed form. One example is the period- 
ically driven two-level system with B{t) - (2y cos 2ojt, 2V sin 2cjt, 2A), or 



//(f) = 



A Ve-^'"^' 



(49) 



where A, V, w e K. The exact propagator is given by 
U(t,0) 



A - cj 

e""^'(cos Of - i ^ sin Qf) 



Q. 

V ■ 
i—e'"' sin Qf 
Q 



V ■ \ 
-i— e"'"' sin Qf 

A - LO 

e"^'(cos Qf + i ^ sin Qf) 



(50) 



with Q - ^J(A - lj)^ + V^. We note that, in accordance with Floquet theory for periodically driven systems, U(7m /a),0) ■ 
Uijijco, 0)" for integer n. The transition probability spin up «-> spin down 



F(f) = |f/2i(f,0)|2 = (^) sin^Qf 



(51) 
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Figure 2: Effective error constant c of 4th-order CFETs with 3 exponentials, given as a function of tlie free parameter /2,i in Eq. I4H . As explained 
in the text, the driven two-level system (Eq. 1491 ) is propagated over the time < t < T = IQtt/lo, for parameters to = I, A = 0.5 and V = 0.5, 1.0. 
Shown are results for CFETs including the A3 term for the optimal choice /2,3 = -0.28 (dashed red curve), and without the A3 term (/23 = 0, solid 
black curve). The horizontal dashed gray line gives c for the CFET CF4:2 with 2 exponentials (Eq. I |38I ). corresponding to 724 = /23 = 0. 

is typical for a Breit-Wigner resonance. 

In the case of two-level systems, application of a CFET requires evaluation of matrix exponentials e"', which 
correspond to propagation with fictitious constant magnetic fields. Each exponential can be evaluated in closed form 
with the relation 

(p _, _, <p 4' -, 
exp(i— n ■ o") = cos — -H i sin —n ■ cr (\n\ - I) (52) 

for the spin 1/2 rotation operator. 

To quantify the CFET error we calculate the deviation e(t) - \\U(t) - U{t)\\ of the approximate propagator U(t) 
from the exact U(t). We use the Frobenius norm for a L x L square matrix 

\\U - Of = ^tv[(U - U)HU -m^^Yj l^'V - U,j\^ , (53) 

ij 

where tr[ ] denotes the trace. This choice is particularly convenient for the Schrodinger equation, where the propaga- 
tors are unitary such that \\U\\ - 1 and \\U - Uf = 2||1 - U W\\. Notice that this definition accounts for phase shps of 
the propagators. With the BCH formula we find that the CFET error 

e{6t) = \\U{St) - U{6t)\\ = 2||1 - e-"<*>e"<*'|| = 2||1 - = 2\\x\\ + 0{5t^^^) (54) 

is indeed determined by the error term x - ^{St) - Q.(6t). The commutator relations of the spin algebra imply that 
the nested commutators Ck in Eq. (l40l i are not independent. This allows us to check the theoretical error analysis 
from Sec. |5] in a situation where cancellation between different Ck plays a role. We note that cancellation is not 
a peculiar consequence of the small Hilbert space of the present example, but of commutator relations dictated by 
physics. Similar relations hold in any relevant situation. 

6.3. Fourth-order CFETs 

To determine the effective error constant c, we propagate the driven two-level system Eq. (|49^ over 20 periods of 
the driving field, i.e. up to a time T = lOnjoj. From the maximal propagation error e = max{e(f)|0 < f < T) we get 
the effective error constant as c = (s/6t)(e/Ty^'^ in the limit 6t « T. 

In Fig. |2] we show c as a function of the free parameter /2,i used in Sec. l5.2l for optimization of 4th-order CFETs 
with 3 exponentials. In both cases (upper and lower panel) c is minimal for /2.1 - 0.45 - 9/20, which confirms 
the previous theoretical analysis based on Fig. [1] In comparison to CF4:2 with 2 exponentials, which has larger c 
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Figure 3: Effective error constant c for unoptimized (CF4:2, CF4:3 from Table|2) and optimized (CF4:30pt from Table|5) 4tli-order CFETs, as 
indicated. As for tlie previous figure, tlie driven two-level system is propagated over T = lOrr/a/. Resuits are given as a function of A, witli dj = I 
and V = 0.1, 0.5, 1.0, 2.0 as indicated in the panels. 



than CF4:30pt, we see that the error reduction is sufficiently large to outweigh the increased eff'ort arising with an 
additional third exponential. We conclude that the optimization is successful and results in more eflicient CFETs. 

The importance of including the A3 term becomes evident when dropping it, i.e. setting /2,3 = (solid black 
curve). Generally, c for such CFETs is large because of significant contributions from the A3 terms in Eq. ( l42l i. and 
the 'optimal' value f2j = 9/20 does not reduce the error Accidental cancellation of different terms occurs for certain 
parameter combinations and leads to the 'dip' in c for /2_i ~ -0.4 (lower panel). Notice that in contrast to such 
artificial minima the true optimized value f2^i ^ 0.45 gives a stable minimum of c. 

In Fig. [3] we show c for a range of parameter combinations of the driven two-level system. Again we see that the 
optimization of CF4:30pt is successful and results in smaller values of c. As an estimate, CF4:30pt is about 10% to 
50% more efficient than CF4:2. Optimization attempts without the A3 terms (CF4:3) result in reduced efficiency. 

6.4. Sixth-order CFETs 

In Fig.|4]we show the effective error constant c for different 6th-order CFETs. We can draw similar conclusions 
as for the 4th-order CFETs. Since the parameter /i 1 = 0.2 of the CFET CF6:5b from Ref. [8] is close to the optimal 
choice /i.i = 0. 16 of CF6:5, both CFETs are comparable, with a slight advantage for CF6:5. The CFET CF6:4 is much 
less efficient, although it requires only 4 exponentials. The optimized CFET CF6:50pt is generally the most efficient, 
while dropping the A4 term (as in CF6:5, CF6:5b) reduces the efficiency. Notice that CF6:5Imp, including the A4 term 
into CF6:5, is not as efficient as the fully optimized CF6:50pt, but still significantly better than the other CFETs. The 
additional freedom of choice of parameters for 6 exponentials (CF6;60pt) does not lead to further reduction of c. 

6.5. Comparison of CFETs of different order 

According to Eq. (|46] |. time-propagation with smaller error, i.e. higher accuracy demands, is more efficient using 
higher-order CFETs. A given Mh-order CFET is most efficient in a certain 'accuracy window', whose size depends 
on the respective error constant c and propagation time T. The intended accuracy goal thus suggests a preferential 
choice of and the corresponding optimized CFET. 

Consider two CFETs of order A^i < N2, with effective error constants ci, C2- Inverting the effort-error relation 
from Sec. 16. II we find that the A^i -order CFET is more efficient than the A^2-order CFET if 

^C2' T 

The decisive quantity is the ratio e/T of the maximal acceptable error e and the propagation time T. 

For a rough estimate, let us assume that the effective error constants ci, C2 are given by the number si, S2 of 
exponentials. With i = 1 , 2, 5, 11 for = 2, 4, 6, 8 we find the following values, which provide some orientation: 
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Figure 4: Effective error constant c for different unoptimized (CF6:4, CF6:5, CF6:5b from Table|3} and optimized (CF6:5Imp, CF6:50pt, CF6:60pt 
from Table |6) 6th-order CFETs, as indicated. The solid black (red) curve corresponds to CF6:5 (CF6:50pt). The propagation parameters are 
identical to Fig. [5] 



en-or elT: ... 6 x 10'^ ... 2 x 10'^ ... 6 x 10''' . . . 
favourable A^: 2 | 4 | 6 | 8 

As a rule of thumb, the accuracy window spans three orders of magnitude: 4th-order CFETs are good for low (error 
10"^), 6th-order for moderate (error 10 "^), and 8th-order for high (error 1 accuracy demands. The use of 2nd-order 
CFETs such as the midpoint rule should be avoided. Long propagation times shift the advantage towards higher-order 
CFETs. 

For a case study we show in Fig. |5]the error-effort plot for 2nd- to 8th-order CFETs, applied to the two-level 
system from Sec. 16.21 for short (left panel) and long (right panel) propagation time. Notice that a very small error 
can be achieved before it saturates at machine precision. The accuracy window of the A^th-order CFET is bounded by 
the crossing with the N + 2 curves. For a moderate error 10"^ (the square root of machine precision for FORTRAN 
double precision numbers), switching from the 4th- to the 6th-order CFET reduces the effort by a factor of 2-3. For 
longer propagation times (right panel) the accuracy window shifts to larger errors, and the 8th-order CFET becomes 
more efficient. The performance of the 2nd-order midpoint rule is several orders of magnitude worse. To illustrate the 
benefit of optimization, we include results for the unoptimized CFET CF6:5. As can be seen, it is never competitive 
in comparison to the (optimized) 4th- or (unoptimized) 8th-order CFET. 

6.6. Time-step selection 

In practice a prescribed accuracy goal has to be achieved without knowledge of the exact solution of the problem. 
A simple, conservative approach is to perform calculations with an ever decreasing time step 6t until convergence, i.e. 
two subsequent calculations agree within numerical round-off errors. This approach wastes much computational time 
if we seek less accurate results, as it tries to construct the (numerically) exact solution. 

For a better time-step selection we can use the known scaling of the error as e = cTSt^ . An estimate of the error 
constant c is obtained from two calculations with different time steps 6t[, 6t2 according to the relation 

max ||xi(f) - X2(0II = cT\6t^ - 6t^\ . (56) 

It involves only the difference between the two approximate solutions xi{t), X2(t), but not the unknown exact solution. 
A reasonable choice is St\l6t2 = (2 . . .3)'^'^, such that the error decreases by a small but significant amount. From 
the estimate of c we can extrapolate to the required time step St for the given accuracy goal. For a reliable estimate 
of the final error it is recommended to perform an additional calculation with smaller 6t. An alternative is to compare 
numerical solutions obtained with two CFETs of different order. 

For applications where the time-dependence of A{t) does not change significantly with f, the required time-step 
can be determined for some finite period < f <K T that is characteristic for the dynamical evolution of the system. 
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Figure 5: Error-effort plot for the Mh-order CFETs CF2:1, CF4:30pt. CF6:50pt and CF8: 11 as indicated, applied to the driven two-level system 
Eq. )49t with a) = I, A = 2, V = 0.5. The system is propagated over T = Snjo) (left panel) or T = lOOnjii) (right panel). Also included is the error 
of the unoptimized CFET CF6:5 (dashed green line). 
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Table 7: Points and weights for Gauss-Legendre quadrature over [0, 1] up to order 8, see Eq. | |571 . Generally, xu+\-m = 1 ~ x,,, and WM+\-m = 



The solution over the entire propagation time [0, T] is then computed with the fixed, predetermined value of 5t. In 
other situations, we can proceed similar to heuristic strategies for general differential equation solvers ifTilfpTll . which 
achieve the global accuracy goal through control of the local time-stepping error If the above extrapolation for 6t is 
performed at every step, it allows for propagation with adaptive time-step selection. 



7. Gauss-Legendre quadrature 

In numerical applications the terms A„ from the Legendre expansion Eq. (l25t can be calculated with a numerical 
quadrature formula. For an optimized A^th-order CFET the quadrature formula must be of order -i- 1 . A convenient 
choice is Gauss-Legendre quadrature (v^ with Njl + 1 quadrature points. 

Gauss-Legendre quadrature is specified through M points x\, . . . , xm, which are the zeros of the Legendre poly- 
nomial Pm{x), and weights wi, . . . , wm (see Table|7]l. The integral of a function f{x) is approximated as 

„i M 

I f{x)dx X V w,„/(x„,) . (57) 

Using the orthogonality of Legendre polynomials it can be shown that Gauss-Legendre quadrature is of order 2M, in 
the sense that this expression is exact for polynomials with maximal degree 2M - 1 . Equivalently, the error of the 
approximation ^ f{t)dt a 5f 2m=i w„,f(x,„6t) scales as df^^^. 
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For the integrals in Eq. (IZSl l Gauss-Legendre quadrature with M - N/2 + I points gives 



N/2+l 



An ^ (2n-l)St ^ w,„P„^i(x,„)A(x,„6t) (58) 

m=I 

for the terms Ai, . . . , A^f/i+i of an optimized Mh-order CFET. This expression can be inserted into Eq. (|29l to obtain 

N/2+l 

Q.i = 5f ^ gi,mAix„,6t) (59) 

m=l 

as a hnear combination of A(f) at different times x„,6t in [0, 5f], where the new coefficients are 

N/2+l 

gi,m = vv„, ^ (2n - l)P„-i(x„)fi^„ . (60) 

We note that, using Legendre polynomials, the calculation of the „, from the tabulated fi „ is much simpler than for 
an expansion in powers of 6t (cf. Refs. flSl)- Specifically for the CFET CF4:2 from Eq. (l38T l we obtain 



r~rw r. f3-2V3,(i) 3 + 2V3 .2)li r. 



f3 + 2V3 3-2V3 p)l 



(61) 



where A*" = A[(l/2- V3/6)5f], A^^) =A[(l/2+ V3/6)5f]. 

It remains to show that Gauss-Legendre quadrature with N/2 + 1 points correctly reproduces the Q,-. If we insert 
the expansion Eq. (l24l l into Eq. (l58T l. we find that the A„ are approximated as 



^72+1 



A„ ^ (2n - 1)^A; ^ WmPn-l(Xm)Pl-l(Xm) ■ (62) 



/> 1 m= 1 

The summands on the right hand side are the + 2-order Gauss-Legendre approximations 

>1 N/2+l 



I P„-l(x)Pl-l{x)dx ^ / W,nP„-l{x,„)Pl-l{Xm) 

Jo t~\ 



(63) 

m=l 

of the scalar product of Legendre polynomials. As long as(n-l)H-(/-l)<A^H-l, i.e. n + l < N + 3, the approximation 
is exact and gives the correct value 6„i/(2n - 1). In particular for n = 1, 2, all integrals for 1 < / < H- 1 are evaluated 
correctly, and Gauss-Legendre quadrature constructs the terms Ai, A2 with an error of order 6t^^^, as required for an 
optimized Mh-order CFET. For « > 3, the integrals with / > - (n - 3) are not evaluated correctly, and introduce an 
error of order dt'^^'^^" into the term A„. 

To understand why the CFET order is nevertheless preserved we must revisit the property (P2) discussed in 
Sees. 13.31 14.11 It states that every nested commutator [A„, , . . . , A„„J contributing to Q.(t) fulfills the condition iik < 
1 + Yjii^k for all = 1, . . . , m. By rule (R2) for the CFET construction this property carries over to the approximate 
Q(f) from Eq. OTb . Since the error of the term A„j incurred from numerical quadrature is of the order 5?'^+'*""', the 
error of the nested commutator is of order (5f^+"^""'+2,»i n, (jQg iq ^jjg multiplication with the remaining terms. By the 
above condition this is at least of order 6t^^^, as required. 

We note that the above argumentation shows the intrinsic connection between Gauss-Legendre quadrature and 
the property (P2) about the absence of certain nested commutators from 0(f). The connection is established through 
expansions in orthogonal Legendre polynomials. Of practical interest is that Gauss-Legendre quadrature with N/2 + 1 
points suffices for (optimized) A'th-order CFETs, although in principle most terms A„ are reproduced with an error of 
lower order. 
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8. Implementation 



Due to the simple product form of Eq. (l28T l the appUcation of CFETs is straightforward. The single difficult 
numerical part is the evaluation of the matrix exponentials e"', which is possible with the Krylov technique discussed 
below. Using Gauss-Legendre quadrature each Q,- is obtained from A(f) as a weighted sum (Eq. (|59]l). For large-scale 
problems, where A(t) is a sparse matrix, it implies that also the Q,- are sparse. Moreover the sparsity pattern of A(t), 
i.e. the distribution of nonzero entries, is preserved: Zeros add up to zeros. This allows for seamless integration 
of CFETs into existing programs, which implement specific data storage formats or matrix-vector multiplication 
routines [18]. The extension to time-dependent Hamilton operators requires only minor modifications. The feature of 
easy implementation gives CFETs additional advantage over the original Magnus expansion. 

8.1. Calculation of exponentials 

Two powerful approaches for the computation of matrix exponentials, particularly of e with sparse hermitian 
matrices M, are the Krylov ifloiEoll and the Chebyshev technique (tA^. Both techniques calculate e^^'^ip, the expo- 
nential applied to a vector, iteratively. They avoid diagonalization of the matrix M, which enters only through matrix- 
vector multiplication as required for sparse matrices. If M is a sum of simple terms, split-operator techniques [15] can 
reduce the computational effort considerably. Other methods, such as the 2nd-order Crank-Nicholson approximation 
g-iM - (1 _ iM)/(l + \M), are neither very accurate, nor suitable for large-scale problems ||22|]. 

The Chebyshev technique is based on the expansion of the exponential function e"' in a series of Chebyshev 
polynomials. Similar to the calculation of spectral functions |23], it has the advantage of low memory demands, 
simple implementation, and unconditional stability and concomitant accuracy for arbitrary large propagation times. 
The main disadvantage, especially for time-dependent Hamilton operators, is the need to determine a-priori bounds 
on the eigenvalues of the matrix M. 

The Krylov technique is based on the Lanczos iteration. Starting with the initial vector i/r, each multiplication with 
M gives a new vector M^\jj, which is orthogonalized to the previous vectors from the iteration. A few K iterations 
generate an orthogonal basis of the Krylov subspace spanned by the vectors ij/. Miff, M^ifr, . . . , M^^^tjj. The exponential 
e'^^ij/ is approximately evaluated within the low-dimensional Krylov subspace, which effectively reduces the problem 
to the calculation of an exponential of a dense m x m matrix 1I22I1 . The success of this procedure depends on the quality 
of the Krylov approximation of M. For the calculation of the exponential e"*^, the error bound 



(ff ^2p</:), (64) 



error < const, x e '''^'^ 

where the constant is independent of K and p, can be established 12(J. Here, 4p = /l^ax - ^^min is the spread of 
the maximal and minimal eigenvalue /Imax, '^mm of M. Increasing K leads to fast reduction of the error However, 
the finite main storage restricts the size of K. Therefore, the Krylov technique requires time-stepping, based on the 
equality g^'^"*' = (e"'^*')", if the error is not sufficiently small after a single Lanczos iteration. For fixed K, the Krylov 
approximation error for g^'™' is of order 6t^ . 

8.2. Comparison of the Krylov and Chebyshev technique 

In Fig.|6]we compare the Krylov and Chebyshev technique with an mth-order Taylor expansion of the exponential 
e^^'^t//, where H is the diagonal 50x50 matrix with entries H„„ = ncj and the vector elements i(/„ are chosen at random 
(prior to normalization of ip). This corresponds to time-propagation for the quantum harmonic oscillator (cf. Sec.|9j. 
The error is given by the ^^-norm e = |i/'(f) - (^(f)! between the numerical result ij/(t) and the exact i//{t), here with 
tlr„{t) - e^"'^"i//„. The effort is equal to the number of evaluations of Hifr in the computation. For the Chebyshev 
technique and Taylor expansion, which evaluate the exponential at once, this corresponds to the number of terms kept 
in the respective series. These definitions of error and effort are also used in the following examples. 

For the left panel in Fig. |6] the system is propagated for f - n/ilQai), i.e. a 20th of the oscillator period. The 
plot shows the typical problems of the Taylor expansion, whose error grows initially before it saturates far above 
machine precision. Since the Taylor expansion violates unitarity, the large errors of the exponential spoil the stability 
of time-propagation. Notice that an Mh-order Runge-Kutta method applied to e is equivalent to using the Taylor 
expansion, which explains their diminished usefulness for quantum systems. For the Chebyshev technique, the error 
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Figure 6: Comparison of the Krylov and Chebyshev technique with the Taylor expansion for the calculation of the exponential e^^^ijj, for the 
Hamilton operator of the harmonic oscillator (//„„ = ntS) as explained in the text. Left panel: Propagation over a 20th oscillator period (; = 
^■/(lOixi)). Right panel: Propagation over a 100th {t = ff/(50(ii)) and a full oscillator period (; = Inlcn). Notice the asymptotic decay of the Krylov 
error ~ effort'^"'. 



decays fast after the first 10 - 20 terms. Unitarity is again achieved only at the level of machine precision, which 
however now can be reached easily. The Krylov technique is competitive for sufficiently many Krylov vectors in the 
iteration {K > 10) and moderate accuracy demands. Notice that the eigenvalues of the quantum harmonic oscillator 
occur as multiples of co, which leads to a large eigenvalue spread p in Eq. (|64] |) and increases the computational eff'ort 
more than the 'classical' time-scale l/co may suggest. 

The Krylov technique becomes more efficient for small time-steps St. It complements the Chebyshev technique 
which excels for longer propagation times. Both scenarios are depicted in the right panel in Fig. |6] This makes the 
Krylov technique the more suitable choice for combination with CFETs, where the length of the time-step is restricted 
by the time-dependence of A{t) (or H(t)). Its central advantage, however, is that it strictly preserves unitarity even 
for finite error This allows us to dispense with the evaluation of the exponentials e^' to very high accuracy when the 
overall error is dominated by the CFET error. Instead, we can use the Lanczos iteration with small K (it must K > N 
for an Mh-order CFET). The reduction of the time-step 6t, in order to decrease the CFET error, reduces the Krylov 
error at the same time. While we recommend the use of the Krylov technique for CFETs we must also note that the 
present example shows that the Chebyshev technique should not be finally dismissed even for short-time propagation. 



9. Example: CFETs applied to the parametric harmonic oscillator 

A genuine example for driven system is the quantum parametric harmonic oscillator 

Hit) = \f + , (65) 

where we allow for a time-dependent oscillator frequency (y(f) - o)^^ + ^ cos Qf . Position q and momentum operator 
p obey the canonical commutation relation [q,p] - i. The oscillator position q(t) - {q)(t) = {tl/{t)\q\i//(t)), given as the 
expectation value of q, follows the classical equation of motion - the Mathieu equation - 

q + (ajl+ ^cos Q.t)q = . (66) 

9.1. Classical oscillator 

The solution of the Mathieu equation provides us with the classical propagator U(t, 0), which is a 2 x 2 matrix. 
According to Floquet theory, the eigenvalues A of C/(27r/0), the propagator over one period, determine the stability of 
the classical system: It is stable, i.e the solutions of Eq. ( |66] | are bounded, if all \A\ < 1, and unstable otherwise. The 
left panel of Fig.|7]shows the stability chart of the parametric oscillator, which we obtained with the CFET CF6:50pt. 
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Figure 7: Left panel: Stability chart of the parametric harmonic oscillator, and the associated Mathieu equation Eq. {66). In the shaded regions 
of parameter space exponentially growing solutions exist, i.e. the system is unstable. For f — ► 0, isolated instabilities occur at integer values of 
2(t)o/n. Right panel: Error-effort plot for the classical parametric harmonic oscillator, propagating a solution cj(t) of Eq. )66t for (oio/il)^ = 2, 
f/tf = 1 over < ? < r = 20n'/(i), i.e over 10 periods of the frequency modulation. The insets display y(t), with initial condition q(0) = 3, q = 0, 
for the same f/n^ = 1 and {ojo/af = 2 (stable regime, left inset) or (pJo/il)^ = 1 (unstable regime, right inset). 

In the right panel of Fig. Qwe show the corresponding error-effort plot for CFETs of different order, where the 
error e = max \q{t) - q{t)\ is measured as the difference between the exact and numerical position q{t) and q(t) over 
10 periods < f < r = lOn/Q.. The optimized 6th-order CFET CF6:50pt is advantageous for practical accuracy 
demands. In the left panel of Fig. [8] we compare different CFETs over a range of ^, values. Shown is the effort 
needed to achieve an error of 10 *' or better As expected, the CFET CF6;50pt is the most efficient. 

9.2. Quantum oscillator 

For the quantum oscillator, we represent position q - {b + Z? ' )/ V2wo and momentum operator p = i VdLto/2(fe ' - b) 
through bosonic ladder operators [b, b''] = 1 . The Hamilton operator is given by 



For ^ = 0, with = wo, we recover the standard Hamilton operator H - ojo(b*b + 1 /2). Truncation of the infinite- 
dimensional bosonic Hilbert space, excluding high energy states, is required to obtain the Hamilton operator as a 
matrix. For the examples we keep the lowest A^;, = 50 Fock states |n), with b'^b\n) - n\n). 

In the right panel of Fig. [8] we compare different CFETs for the quantum oscillator As the initial wave function 
we choose a coherent state \ijjc) with (i/'c|§|i/'c) = 3, {ijjAPVI^c) - 0. The error is measured by the deviation e - 
max |i/r(f) - ^(f)| between the exact and numerical wave function i/'(f) and (^(f), over 10 periods < f < 207r/O. Shown 
is the effort needed to achieve an error of 10"^ or better, as for the classical case. 

A new aspect in comparison with the CFET error analysis for the classical oscillator is the numerical evaluation 
of the exponentials e"' with the Krylov technique. For few Krylov vectors {K - 10, upper panel) the Krylov error 
from the approximate exponentials dominates. In this case, short time-steps are preferential to reduce the Krylov 
error sufficiently, with the consequence that the unoptimized 4th-order CFET CF4:2 is most efficient since it uses the 
smallest number of exponentials. With more Krylov vectors {K - 20, lower panel) the exponentials are evaluated 
to much higher accuracy also for longer time steps, and the expected advantage of optimized higher-order CFETs is 
recovered. The overall most efficient propagation is achieved with the CFET CF6:50pt for K = 20. 

Notice that the necessary Hilbert space truncation limits calculations in the unstable regimes shown in Fig.|7] as 
the classical instability manifest itself for the quantum system in the excitation of high energy Fock states. Although 
the truncated Hamilton operator remains hermitian and can be used for time-propagation, the position expectation 
value {q){t) cannot be expected to obey the classical equation Eq. ( |66] |. 




l)(/,t2 + /,2) + (f^ + l)(2/,l7,+ l) 



(67) 
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Figure 8: Effort required to acliieve an error 10"^ for time-propagation of the classical (left panels) and quantum (right panels) parametric harmonic 
oscillator Eq. | |651 over < ? < T = lOnja), i.e. 10 periods of the frequency modulation. Solid (dashed) curves corresponds to optimized 
(unoptimized) CFETs as indicated. For values in the gray shaded regions, unstable solutions of Eq. l l66t exist according to Fig.|2] Left panels: For 
the classical oscillator, the effort is given as the number of evaluated exponentials. Shown are results for (aJolSl)^ = 2 as a function of ^jOr (upper 
panel), and for ^/£1^ = 1 as a function of (oJolQ.)^ (lower panel). Right panels: For the quantum oscillator with ^/£1^ = 1, the initial coherent 
state is propagated using Nj, = 50 Fock states. The effort counts the number of applications of the Hamilton operator The exponentials are 
evaluated with the Krylov technique, for K = \0 (upper panel) and K = 20 (lower panel) Krylov vectors. The green curve in the upper panel 
reproduces the result for CF6:50pt, K = 20 from the lower panel, as indicated. 



9.3. The interaction picture for numerical time-propagation 

Standard time-dependent perturbation theory is based on the interaction picture. Consider a decomposition A{t) - 
D + B{t), where Z) is a constant diagonal matrix. The propagator for D is the exponential e'^. The interaction picture 
is defined by x'{t) = e^"^x(t). If B(t) = 0, x'{t) is constant. Otherwise, it obeys the equation of motion 

d,x'(t) = e-'"B(t)e'"x'(t) = B'(t)x'(t) , (68) 

where B'(t) = e"'^B(f)e'^. Since D is diagonal, the matrix elements of B'{t) are easily calculated, with - 
g(Am-A,,m)'/j/^^ Notice that the diagonal elements of B(t) do not change, and a sparsity pattern is preserved. 

The interaction picture is useful if it simplifies the equation of motion when B{t) is a small perturbation. That is is 
generally not the case can be understood for the driven two-level system Eq. ( l49l l from Sec. 16.21 where D is identified 
with the term Act,. The equation of motion in the interaction picture is identical to the original equation of motion 
with new parameters A' = 0, = w - A. As can be seen from Eq. (ISOl l, the propagator in the interaction picture 
is identical to the original propagator apart from an additional rotating phase e'^"^'. This implies that the calculation 
in the interaction picture has not simplified. From the perspective of numerical time propagation the difficulty even 
increases since B'{t) varies faster than B{t) due to the additional time-dependence acquired in the transformation with 
e'^. This is particularly true if B{t) is a small perturbation, since - A\ > oj for large A. 

Notice that for the present problem the choice D = ±wcr, leads to a constant Hamilton operator in the interaction 
picture, which allows for the construction of the exact propagator Eq. ( l50t . Indeed, the celebrated rotating wave 
approximation is exact for this particular case. Despite its persistence in quantum optics it does not easily generalize 
to other situations. 



9.4. The interaction picture for the harmonic oscillator 

While the interaction picture per se does not simplify time-propagation, it can be useful to reduce the computa- 
tional effort associated with the numerical evaluation of exponentials. As discussed in Sec. 18.21 the quantum harmonic 
oscillator is an example where large eigenvalues nu) increase the effort. Switching to the interaction picture, with 
D - (ob^b, increases the CFET error because of the additional time-dependence on the scale of co, but simplifies 
the evaluation of exponentials since the large diagonal entries no) are removed from the matrix. We illustrate this 
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Figure 9: Time-propagation in the interaction picture (IPIC) of tlie quantum parametric liarmonic oscillator Eq. ( 1651 with ^/£1^ = 1 over < ? < 
T = lOn/oj. Left panel: Error-effort plot for (oi/Cl)^ = 4, using the CFET CF6:50pt and the Krylov technique with K Krylov vectors as indicated. 
The solid curves give the error from standard propagation, the dashed curves from working in the interaction picture. The gray dashed line shows 
the asymptotic scahng of the CFET en'or ~ efforF*. Right panel: Effort required to achieve an en'or 10"* with the CFET CF6:50pt, as a function 
of (wo/n)2. For comparison, results for standard propagation with K = 20 of Krylov vectors are reproduced from the lower right panel in Fig. [8] 



possibility with the error-efFort plot for the CFET CF6:50pt in Fig. |9](left panel), where results for the interaction 
picture are compared to those from standard propagation for a different number K of Krylov vectors. We see that in 
the interaction picture the Krylov error is much reduced so that the CFET error, with scaling ~ effort^, dominates over 
the entire range. For moderate accuracy demands, with errors down to 10"^, the interaction picture with only K - 1 
Krylov vectors is most efficient. For smaller error, the interaction picture is again less favourable, since the CFET 
error has increased in comparison to standard propagation. The right panel of Fig. |9] shows the effijrt to achieve an 
error 10"^, supporting the expectation that the interaction picture becomes rather efficient at larger {(dqIQ)^. Whether 
there is a benefit of using the interaction picture also for non-bosonic systems remains to be studied. 

10. Comparison of CFETs to Floquet approaches 

For problems with a periodic time-dependence Floquet theory suggests exploitation of the periodicity of the prop- 
agator A notable implementation of this idea is the (f, f')-method |24]. Introducing time as an additional variable f', 
the wave function if/if) is recovered from the solution f') = exp(-i'Kf)*P(0, f') of the Schrodinger equation with a 
time-independent Hamilton operator - H(t') - idf as i/'(f) - *F(f, f). The validity of this procedure can be checked 
by evaluation of 5,*F(f, t), with initial condition *F(0, f') = '/'(O). In computations, the auxiliary degree of freedom t' 
is represented with a Fourier basis of periodic functions 0„(f') - g^mn/'/r rpj^g calculation of the matrix exponential 
exp(-i'7/f) in the enlarged Hilbert space is ideally suited for the Chebyshev technique providing solutions for one or 
more periods at once. The accuracy is determined by the number Nf of Fourier modes kept in the calculation. 

In Ref. ll24ll . the (f, f')-method was compared to a 2nd-order Magnus propagator It was found that the {t,t')- 
method is far more efficient and allows for reduction of the error down to machine precision with moderate effi3rt. 
Following these examinations, we consider the quantum harmonic oscillator H - 12 + 12 + f(t)q with a time- 
dependent periodic force /(f) = /(f -I- T). We propagate the initial coherent state li/^c) over 10 periods {T - 5/37r) with 
(i) a sinusoidal force /(f) = sin^27rf/r, (ii) a Gaussian pulse /(f) = exp(-((f - r/2)/0.4)^). The error-effijrt plot in 
Fig. [TO]compares the (f, f')-method with higher-order CFETs. 

We see that in both examples the (f, f')-method is significantly less efficient than any but the 2nd-order CFET. 
Although the f-f' error drops rapidly once Nf is sufficiently large to represent the Fourier components of the aux- 
iliary wave function *I'(f, f'), even moderate accuracy requires Nf > 2^ and proportionately large effi^rt. For the 
Gaussian pulse more Fourier modes must be kept, since weight is distributed to higher Fourier coefficients /„ — 
HIT) ^ f(t)e^™'l'^dt of the dr iving force /(f). This restricts the use of the (f, f')-method if memory limitations are 
a concern. Notice that splitting the periodic problem into several time-steps increases the effijrt further, in particular 
since the Fourier coefficients of the then discontinuous force decay more slowly. 
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Figure 10: Comparison of the (t, f')-metliod with CFETs for the forced quantum harmonic oscillator Similar to Figs.[8]|9] the initial coherent 
state \if/c) is propagated in a truncated Hilbert space with Nt = 50 Fock states over 10 periods of the sinusoidal force (left panel) and the Gaussian 
pulse (right panel) specified in the text. The black circles coiTespond to the (t, f')-method with Nf Fouiier modes as indicated, the other curves to 
CFET/Krylov propagation with A" = 10 Krylov vectors. The insets display the Fourier coefficients |/„p and f{t) itself. 



The poor efficiency of the (f, f')-method in comparison to the higher-order CFETs is not a failure of the Floquet 
approach. If we associate a fictitious time-step T/Nf with the representation of the wave function *i'(f, f') through A^^ 
Fourier nodes per period, it is much larger than the time-step in the CFET time-stepping. This is in accordance with 
the expectation that for periodic problems Fourier decomposition provides a better representation of the propagator 
than the concatenation of step-wise constant propagators. A related observation is the increased accuracy of the 
Fourier transform for integration of periodic functions over the combination of finite order polynomial integration 
formulae. In total, the (f, f')-method requires less application of the Hamilton operator 'H for propagation over the 
entire 10 periods than the CFET/Krylov technique with short time-steps. However, the practically relevant effort of 
computations in the Fourier space is just larger by Nf, which, effectively, renders the (f, r')-method less efficient than 
higher-order CFETs. 



11. Further applications 

We complete our study of the practical applicability of CFETs with calculations for two complex quantum systems, 
for which neither exact solutions nor classical analogues are known: A chain of interacting spins - or two-level atoms 
- in pulsed magnetic fields (Sec. lll.fT i. and the hydrogen atom in an electric field (Sec. I11.2l i. Both systems feature 
non-trivial physical effects, and require computation of exponentials for moderate -to-large sparse matrices. 

11.1. Driven spin chain 

In first approximation atoms in a strong light field can be described by interacting spins 1 /2 in a magnetic field. 
We consider the Hamilton operator 

H = j] + / " + (69) 

.5=1 .S=l 

of a spin chain with S spins, where 

//(^^ = AcrW + %V(t) crW - 3 V(t) cr« = (^f^^^ ^^^^ J (70) 

is the Hamilton operator of a single spin, subjected to a magnetic field similar to Eq. (l49T l. If the system is initially 
prepared in the ground state, the magnetic field induces transitions to excited states. For the choice 

V(t) = — — , (71) 
cosh t/T 
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Figure 11: Time-propagation of the driven spin chain Eqs. | |69K I70K with A = 1, / = 0.1. The system is pumped by a sequence of resonant pulses 
(Eq. 17 U with a) = A, t = \,V = 1/(4t)), centered at multiples of fo = 9;r/2. The initial state at -fo/2 is the product state |i . . . I). Left panels: 
Total spin z-component a'-(1) as a function of time, for S = 20 spins. For comparison, the curve for a single spin (S = 1) is included. The upper 
panel shows the pulsed field 51 V(f). Right panel: Error-effort plot for the time-propagation shown in the left panel, for different CFETs and number 
of Kiylov vectors K as indicated. For comparison, the panels include either the curve for CF6:50pt, K = 10, or for CF6:5 in the upper right panel. 



a magnetic pulse of half-width ^ 1 .32t and frequ ency w, the transition probability for a single spin (5 = 1) can be de- 
duced from the result for the Rosen-Zener model |j25|. Specifically, the transition probability Poo - |(t| t/(oo, -oo)|J,)p, 
i.e. the probability that the spin is flipped through the pulse, is 



sin^ jiVt 



cosh^ ;7r(A - (x>)t 



(72) 



In Fig.[TT|(left panel) we show the expectation value a^it) = I S)Yf^^\{4'it)\crf\4>{t)) for a sequence of magnetic 
field pulses. The pulse sequence brings a single spin (curve for S - Y) back to its initial state after two subsequent 
pulses. For several interacting spins (S - 20), dephasing leads to a state with iT;(f) = after the first few pulses. 

The right panel in Fig. [TT]compares the efficiency of different CFETs with a different number K of Krylov vectors. 
This example shows, similar as for the harmonic oscillator, the importance of balancing the Krylov and CFET error 
For small K - 1 the Krylov error dominates, which gives the 4th-order CFET CF4:30pt an advantage over higher- 
order CFETs because it requires less exponentials per time- step. The Krylov error is however less dominant than for 
the harmonic oscillator, and the 6th-order CFET CF6:50pt with K - \Q results in the most efficient time-propagation. 
Notice that the unoptimized CFET CF6:5 (upper right panel) is about 50% less efficient. As an interesting feature we 
note that the slope of the curves for K - \Q resembles that of a 9th-order relation (error ~ 1 /effort''), which is the 
expected scaling of the Krylov error for K - \Q (cf. Eq. (|64]|). The 'bend' from the 9th-order scaling to a 4th-order 
scaling is clearly seen in the curve for CF4:30pt. The error of higher-order CFETs remains smaller than the Krylov 
error, and 9th-order scaling persists down to machine precision. 



11.2. The hydrogen atom in an electric field 

Our last example is that of a hydrogen-like atom in a classical monochromatic electric field along the z-axis. The 
Hamilton operator in dipole approximation is // = - -- + E,(t)d;, where E,(t) denotes the field strength and d, = z 
is the z-component of the dipole operator Working in the basis of hydrogen eigenstates \nlm}, with energy w„ - -l/n^ 
for ^^(f) = 0, the quantum number m is conserved for the above Hamiltonian. We consider only the m - sector The 
required matrix elements of the dipole operator d^ can be calculated analytically or with a one-dimensional numerical 
integration. They are non-zero only between states for which the respective I differs by ±1. 

The system is initially prepared in the E~{t) = ground state i/'(f = 0) = |10), and 4r(t) is calculated for < f < 
T - 10^ using the CFET CF6:50pt in combination with the Krylov technique (K = 10). The electric field is given 
by Er(t) - E''^h{t) cos Q.t, where h{t) = (1 + «)/(! + aexp(-b(x - fo)^)) is an envelope function with a - b - 10"^, 
to - 5000. In Fig. [12] we show the summed occupation probability P„(t) - Y/i=o K"'l'/'(0)P (left panel) and its time 
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t t a 

Figure 12: Time evolution of a hydrogen atom in an alternating electric field. Left panel: Occupation probability P„(t) as a function of time, for 
n = 0.27, £? = 0.1. The upper row gives £,(?)■ On the left we show the envelopes, while the magnification on the right resolves the fast oscillations 
(the curve for n > 4 is omitted). The shown propagation time covers 430 field oscillations. All n < 50 states are kept in the calculation. Right 
panel: Time-averaged occupation probability P„ as a function of f2, for e9 = 0. 1 . The lower panel shows the result for the three level system 1 10), 
|20), |21). The vertical gray lines indicate the resonance frequencies 3 /{4k), k = 1,2,..., for transitions |10) m |21). 

average P„ - (l/T) J^^ P„{t)dt (right panel). In the weak coupling limit E, — > 0, resonances occur if the transition 
frequency |(y„, -<y„,| between states |«i, /) and \n2, l±\) is a multiple of the field frequency O. This behaviour is clearly 
seen if only the three « = 1, 2 states |10), |20), |21) are included in the calculation (lower right panel in Fig.fTSTi. The 
broad resonance at Q = wi - W2 = 3/4 is most pronounced, while the resonances at Q - 3/8,3/12,... become 
increasingly sharp (for a non-classical field, these would correspond to multi-photon absorption). Inclusion of states 
with larger n (upper right panel, with n < 50 in the numerical calculation) shifts the frequencies of the n = 1 <-> n = 2 
transition, and leads to the numerous sharp resonances of transitions to higher excited states. 

12. Conclusions 

The development of practicable techniques for the propagation of driven quantum systems requires realization of 
high theoretical efficiency gains under the restrictions of actual applications. In the present paper we studied a par- 
ticular class of numerical techniques, the commutator-free exponential time-propagators, which combine favourable 
theoretical properties, such as preservation of unitarity and high approximation order, with the virtue of simple imple- 
mentation. 

Conceptually, CFETs are related to the more traditional Magnus expansion. From the practical point of view, they 
are in fact the better alternative, at least for the problems studied here. Avoiding commutators makes them easier 
to implement and also more efficient, since the complicated structure of the original Magnus expansion and all the 
bookkeeping it requires is replaced by their simple exponential product form. 

We dealt with the derivation, optimization, and application of CFETs from the common point of view of the prac- 
titioner who wants to solve the Schrodinger equation. For every issue the present work extends the existing literature. 
Our construction and analysis of CFETs relies essentially on the use of Legendre polynomials and their orthogonality 
properties. In this way we can provide a comprehensive and self-contained presentation. It also simplifies the error 
analysis and allows us to identify the importance of including higher-order terms for the CFET optimization. We 
provide coefficients of fully optimized 4th- and 6th-order CFETs, as well as of a good albeit unoptimized 8th-order 
CFET. As both the theoretical and practical error analysis show full optimization is successful in further reducing the 
error, leading to about 50% higher efficiency in comparison to the partly optimized counterparts. While the potential 
of 6th-order CFETs is probably largely exhausted, optimization of 8th-order CFETs remains promising. 

We have discussed the practical application of CFETs at great length, paying particular attention to realistic situ- 
ations where exponentials can not be calculated in closed form. Based on our findings, we generally recommend the 
use of the CFET CF6:50pt together with a Krylov calculation of the exponential using about 10-15 Krylov vectors. 
The results for the examples presented show that very accurate results can be obtained with moderate effort. They 
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provide evidence that for the Schrodinger equation optimized higher-order CFETs are substantially more e fficient than 
alternative techniques such as general purpose Runge-Kutta methods or numerical Floquet approaches. Most impor- 
tantly, CFETs are robust: They are unconditionally stable, and their quality does not substantially decline at points of 
resonance. CFETs are thus a good choice for library routines for time-propagation. We believe that the implemen- 
tation and optimization of a general purpose time-propagation routine provides most potential for further significant 
efficiency gains. Irrespective of machine dependent implementation details, this has to include refined strategies for 
the automated choice of the step-size and the number of Krylov vectors, as well as tracking of the accumulated error. 
Even now CFETs are a viable and convenient technique for the time-propagation of driven quantum systems. 
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Appendix A. Recursion for the Magnus expansion 

It is possible to write every Q„(f) from the Magnus expansion as an n-fold time-ordered integral 

' 'i 

Q„(/) = J dh J dt2...J dt„Z„(ti,...,Q = J Z„(l, ...,«) (A.l) 

A„[f|l,...,n] 

of a multivariate function Z„(l, . . .,n) = Zn{t\, . . . , f„). With regard to Eq. (flTt . we have 

Zi(fi) = A(fi) , Z2(fi,f2) = ^[A(f,),A(f2)] , Z3(fi,f2,f3) = ^[A(fi),[A(f2),A(f3)]] + i[[A(f,),A(f2)],A(f3)] . (A.2) 

The integration domain of the time-ordered integral is the set of decreasing n-tuples 

A„[f|/i, . . . , /„] = {(fi, . . . , f„) e R"|f > f,-, >?„>■■■> f,„ > 0) , (A.3) 

where /i ,...,/„ is a permutation of 1 , . . . , n denoting the arrangement of the tuple elements. For example, A2 [f 1 1,2] = 
{(fi, f2) e M^|f > fi > f2 > 0) and A2[f|2, 1] = {(fi, f2) e > f2 > fi > 0}. Every permutation selects one of the n\ 
wedge-shaped subsets of the n-dimensional hypercube [0, f]", which is the disjoint union of all these sets (up to points 
from an n - 1 -dimensional subset, which as a set of measure zero is irrelevant for integration). The time derivative 
Q„(f) is given by the n - 1-fold integral 

Q„(f)= J Z„(f,l,...,«-1). (A.4) 

A„_,[f|l,...,n-1] 

We also note that 

Jdf J /(f',l,...,n) = J f(l,...,n+l) . (A.5) 

A„[r'|l, ...,«] A„+i[r|l,...,;i+l] 

According to Eq. (fTST l. Q.„+\ is given as 



a,+i(f) = y f^n^ y r dt'[...[^„,it'),n„,{t')],...,n„jt')],n,^^,{t')] 

t 

= Z J^^' J [z„„z,„...,z„,„j(.M,..., 



(A.6) 



n) 



niH h«,„+i=n+l "III , 7 
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where we use the notation (n - ni + ■ ■ ■ + n^) 

](!,...,«) = [...[Z„,(l,...,ni),Z„,(ni + + «2)], ■ ■ ■ ,Z„j(n - + !,...,«)] (A.7) 

for the nested commutator in the integrand. The integration domain is a product set 

I„[t\ni,...,nk] = A„,[f|l,...,«i] X A„Jf|l,...,n2] X ■ ■ ■ X A„Jr|l, . . . , n^] (n ^ m + ■ ■ ■ + rit) . (A.8) 



To bring the integrals in Eq. (IA.6I 1 into time-ordered form, the integration domain is spht into disjoint pieces that 
are mapped onto the 'wedge' sets A„ [f | . . . ] through a permutation of the integration variables. For every (f 1 , . . . , f„) e 
I n{t\ni, . . . ,nk\ a unique permutation tt exists that orders the «-tuple such that {t„-i(iy . . . ,t„-i(„)) e A„[f|l, . . . ,«]. The 
admissible permutations are those that respect the order of elements corresponding to each of the A„; [f | . . . ] factors in 
I„{t\n\, . . . ,nk\. These form the set 

!P„[ni, . . . , n^] = {tt is permutation of {1, . . . , n)| n{\) < ■ ■ ■ < n{n\) and n{ni + 1) < ■ ■ ■ < ji{n\ + n2) 

(A.9) 

. . . and n(n - + 1) < ■ ■ ■ < n{n)} , 

where still n - n\ + ■ ■ ■ + nk.lt has n!/(ni ! ■ ■ ■ nj^!) elements. In particular, !P„[1, 1, . . . , 1] is the set of all permutations, 
while !P„[n] contains only the identity. 

The decomposition of X„[f|ni, . . . , n^] into disjoint subsets congruent with A„[f|l, . . . , n] is given by 

I„[t\m,...,nk\^ y A„[f|7r-i(l),...,7r-i(«)] . (A.IO) 

;reP„[;ii, ...,«!.] 

Permutation of the integration variables then gives the identity 

J /(!,...,«)= Yj [ f(Ml),..-,n(n)). (A.ll) 



This identity allows us to express the integrals in Eq. (IA.6b as time-ordered integrals. The final integration over f' 
preserves time-ordering according to Eq. (IA.5) . 

After these preparations we can finally state the recursion 

Zi(f)=A(f), 

" / I \m+l 

Z„+i(0, !,...,«) = ^ J— ^ ^ ^ [Z,„,...,Z„,„^,](0,7r(l),...,;r(n)). (^-12) 

m=l «i,...,«,„+i>I ;re?'„[«i-l,...,n„,+i] 

«1 H hrt,j,+ | =rt+l 

While the first terms Z„ can be obtained by hand, the calculation of higher terms is better left to the computer Consider 
exemplarily the calculation of Z3. The sum over ni, . . . , n,„+i contains 2 + 1 terms for m =1,2. Thus, 

Z3(0,l,2)=i 2 [Zi,Z2](0,;r(l),7r(2))+i ^ [Z2,Zi](0,;r(l),7r(2))- ^ ^ [Zi,Zi,Zi](0,;r(l),7r(2)) 

;re!P2[0,2] ;re?'2[I,I] ;re'p2 [0, 1 , 1 ] 

= i[Zi,Z2](0, 1,2) + i([Z2,Zi](0, 1,2)+ [Z2,Zi](0,2, 1)) - i([Zi,Zi,Zi](0, 1,2) + [Z,,Zi,Zi](0,2, 1)) 
/ 2 O 

= |[0, [1,2]] + i[[0, 1],2] + i[[0,2], 1] - i[[0, 1],2] - i[[0,2], 1] = i[0, [1,2]] + i[[0, 1],2] , 
4 4 4 

(A.13) 

writing [0, [1,2]] = [A(0), [A(1),A(2)]] - [A(fo), [A(fi), A(f2)]] etc. as a short-hand notation. This reproduces the term 
from Eqs. (jAjl. 
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Ai, Az, [Ai,A2], [Ai,[Ai,A2]], [A2, [Ai, A2]], 
[Ai,[Ai,[Ai,A2]]], [A2,[Ai,[Ai,A2]]], [A2, [A2, [Ai, A2]]], 
[Ai, [Ai, [Ai, [Ai, A2]]]], [A2, [Ai, [Ai, [Ai, A2]]]], [A2, [A2, [Ai, [Ai, A2]]]], 
[A2, [A2, [A2, [Ai, A2]]]], [[Ai, A2], [Ai, [Ai, A2]]], [[Ai, A2], [A2, [Ai, A2]]], 
[Ai, [Ai, [Ai, [Ai, [Ai, A2]]]]], [A2, [Ai, [Ai, [Aj, [Ai, A2]]]]], [A2, [A2, [Ai, [Ai, [Ai, A2]]]]], 
[A2, [A2, [A2, [Ai, [Ai, A2]]]]], [A2, [A2, [A2, [A2, [Ai, A2]]]]], [[Ai, A2], [Ai, [Ai, [Ai, A2]]]], 
[[Ai, A2], [A2, [Ai, [Ai, A2]]]], [[Ai, A2], [A2, [A2, [Ai, A2]]]], [[Ai, [Ai, A2]], [A2, [Ai, A2]]] 

Table B.8: The 23 Hall basis elements with generators Ai, A2 and up to 5 commutators. 



Appendix B. Free Lie algebras and Hall bases 

Avoiding formal definitions, the basic concept of a free Lie algebra can be understood in simple terms. For more 



thorough accounts, see Refs. Ill 1111211 . 

A free Lie algebra is a vector space equipped with a function in two arguments [■, ■], the commutator. It consists 
of all nested commutators of the generators A i , A2, . . . and all linear combinations thereof. In addition to the standard 
vector space properties, one demands bilinearity [X + Y,Z] - [X,Z] + [Y,Z], [cX,Y] - c[X,Y] and anti-symmetry 
[X, Y] = -[YX] of the commutator, together with the Jacobi identity [X, [Y,Z]] + [Y, [Z,X]] + [Z, [X, Y]] = 0. No 
further relations hold: Two elements of the free Lie algebra are different if they cannot be transformed into each other 
with these identities. In other words, only the minimal relations characteristic for a commutator hold. 

Anti-symmetry and the Jacobi identity imply linear dependencies between nested commutators of the generators. 
In particular, they do not form a vector space basis of the free Lie algebra. For three elements X,Y,Z exist 12 
commutator combinations 

[X, [YZ]], [X, [Z, Y]], [Y, [X,Z]], [Y [Z,X]], [Z, [X, Y]], [Z, [YX]], 
[[X,Y],Z], [[X,Z],Y], [[YX],Z], [[YZ],X], [[Z,X],Y], [[Z,Y],X] . 

Any three of them are linearly dependent, such that we must select two for a basis. With this in mind, the Hall basis 
construction defines a systematic selection rule. First, define an order "<" on the generators and nested commutators. 
For the generators, set A, < Aj if / < j. For the commutators, set [X, Y] < [V,W] if X < V or X - V,Y < W. Set 
generally X < Y if Y is composed out of more commutators than X. The Hall basis is now defined recursively: (HI) 
All generators A, are in the Hall basis, (H2) a commutator [A,-, A^] is in the Hall basis if A, < Aj (i.e. / < j), (H3) if 
X,Y,Z are in the Hall basis, so is [X, [Y, Z]] provided that [Y, Z] is in the Hall basis and F < X < [Y, Z]. 

To understand rule (H3), observe first that it removes the ambiguity due to anti-symmetry, since it enforces X <Y 
for Hall basis elements \X, Y]. Now consider a nested commutator \X, \Y, Z]] from the Hall basis. It is F < X < [Y, Z] 
by (H3), and also Y < Z. Consequently, Y < [X,Z]. Both properties rule out most commutators from Eq. (IB.l) 
apart from [X, [Y,Z]] itself and [Y, [X,Z]], [Z, [Y,X]], [[Y,X],Z]. If X = Y, only [Y, [X,Z]] is non-zero. Otherwise, for 
Y <X, [Y, [X,Z]] violates (H3). Then, depending on whether Z ^ [YX], either [Z, [Y,X]] or [[F,X],Z] fulfills (H3). 
If Z = [X, Y], both commutators vanish. In any case, at most one commutator from Eq. (IB.l) is a Hall basis element 
in addition to [X, [Y,Z]]. This argument implies linear independence of the basis elements, and can be turned into an 
inductive proof. Moreover, rule (H3) amounts to a recursive algorithm to check for membership of the Hall basis. 

Completeness of the Hall basis can be shown with a similar argumentation. Based on this, a recursive al- 
gorithm can be devised to express commutators [X, Y] as linear combinations of the Hall basis elements. In Ta- 
ble lB.81 we show the first 23 Hall basis elements involving the generators Ai, A2. For example, the last commutator 
[[Ai, [Ai,A2]], [A2, [Ai,A2]]] fulfills (H3) with X = [Ai, [Ai,A2]], F = A2, Z = [Ai,A2]. Another example is to write 
the four-fold nested commutator [Ai, [A2, [Ai, [A2,Ai]]]] = [[Ai,A2], [Ai, [A2,Ai]]] - [A2, [Ai, [Ai, [Ai,A2]]]] as the 
unique sum of two elements from the table. As discussed in Sees. [3] 111 only a small subset of all Hall basis elements 
needs to be considered for the Magnus expansion or CFET construction. 
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Appendix C. Order conditions for 6th-order CFETs 

The order conditions for 6th-order CFETs can be largely solved by algebraic manipulations. For 6th-order CFETs 
with 5 exponentials, one has 7 equations for the 8 coefficients /i,i,/i,2,/i,3,/2,i,/2,2,/2,3,/3,i'/3,3' follows: 

Ai: l=2/u+ 2/2,1 +/3,i 
A3 : = 2/1,3+ 2/2,3 +/3,3 

['4l,'42] : - -/. 1/1,2 - 2/2,1/1,2 - /2,l/2,2 - .h,\f\,2 - /3,l/2,2 

o 

[A2, A3] : -— - f\2f\3 + 2/1,2/2,3 + /l,2/3,3 + 12,2/2,3 + /2,2/3,3 

11 2 1 1 1 1 

[^^1, ['4i,'43]] : 77: - +0/1,1/2,1/1,3 - 0/1,1/2,1/2,3 - 0/1,1/2,1/3,3 + 7/1,1/3,1/1,3 - 0/1,1/3,1/2,3 - 7/1,1/3,1/3,3 
oU J 3 3 O 3 O 

1 2 1 2 2 1 1 2 2 1 2 

~ 3/1,1/2,3 - ^/i, 1/3,3 + -/2,l/3,l/l,3 + -/2,l/3,l/2,3 - -/2,l/3,l/3,3 + 3/2,1/1.3 - ^/2,l/3,3 

1 2 1 2 

+ g/3,1/1,3 + g/3,1/2,3 

[A2, [Ai, A2]] ■ = ~^h,'^f\,2 ~ 1/2,1/1,2/2,2 - l/2,l/i^,2 - 3/2,1/2^,2 ~ 1/3,1/1,2/2,2 - 2/3,1/u ~ 2^^,^f2,2 



+ ^/l, 1/3,1/2,2 + 3/1,1/2,1/1,2 + ^./l, 1/2, 1/2,2 + ^./l,l./3,l./l,2 + ^/l, 1/3,1/2,2 + Y^/l, 1/1,2 

(C.l) 



+ J^f2.lflif2.2 + -/|,i/3,l/2,2 + — /2j/2,2 



The order conditions for 6 exponentials have a similar structure, but are too long to be shown here. 

Apart from degenerate cases, the order conditions can be reduced to a single polynomial equation. We consider 
/i,i as a free parameters. Then, if /2,i is the solution of p(/i,i,/2.i) = with the polynomial 

p(x,y) = - 2 + 30x - 192x2 + 680x^ - 1440x'^ + 1815x^ - 1250x^ + 360x^ 
+ (18 - 232x + 1230x2 - 3440x^ + 5345x'* - 4350x^ + 1440x^)y 

+ (-60 + 650x - 2740x2 + 5655x^ - 5710x'* + 2250x^)/ + (90 - 800x + 2535x2 - 3450x^ + niOx'^)y^ 



+ (-60 + 425x - 920x2 + 630x^)/ + (15 - 80x + 90x2)/ 



of degree 5 in y, the remaining coefficients are given by 



(C.2) 



l+5/i,i(/i,i-l) ^ 1 - 6/2,2 + 12/,i/2,2 + 6/2,1/2,2 

/2,2 - :J7T77 r~? 7T77 — , n r r~? 7T ' /i,2 ~ 



/l,3 



./2,3 



30(/i,i + /2,i - l)(/,i + /2,i)(2/i,i + /2,i - 1) 6(1 - /,i) 
(2/i,i-l)(2/i,i+/2,i-l)- 3/2,2 

30(/i,2(2/i,i - l)(2/i,i + /2,i - 1) + /2,2(1 + 8/2i + 2(/2,i - 2)/2,i + (8/2,1 - 7)/i,i)) 

/l,l + 3/1,2 + 4/.1/2.1 + 2(/2.i - l)/2,i + 6/2,2 - 1 

30(./,2(2/,i - l)(2/,i + /2,i - 1) + /2,2(1 + 8/2i + 2(/2,i - 2)/2,i + (8/2,1 - 7)/i,i)) 



(C.3) 



/3,1 - 1 - 2/1,1 - 2/2,1 , ./3,3 - -2/l,3 - 2/2,3 ■ 

Several solutions exist with simple explicit expressions for the coefficients, such as the ten solutions shown in 
Table IC.9I Unfortunately, none of these is competitive with the CFETs from Tables [3] |6] For the CFET CF6:5 
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6th-order, 5 exponentials 



/i.i ^ (5 - V5)/10 /2,i ^ (23 - 4 V5)/60 



/u 


= (5 + V5)/10 


/2,1 


= (23 +4V5)/60 


/u 


= (65 + Vl005)/90 


/2.1 


= 3/10 




= 3/10 


hi 


= (553 ± 3 V2()T)/2400 




= 1 


/2,1 


= (30+ V290±50V5)/60 



Table C.9: Explicit simple solutions of the order conditions Eq. lC.ll for 6th-order CFETs with 5 exponentials. The remaining coefficients can be 
found with Eq. )C.3) . In the last row, all four combinations of the signs are allowed. 



with fi x - 0.16 = 4/25, the coefficient /2,i — 0.387524052 ... is the single real root of the polynomial p{x) - 
-126131602 + 1646347450X - 7919062500^2 + 1695003 1250x^^ - 15834375000x'* + 5498046875x^ 
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